class Gas(object):

    def __init__(self, start_date: datetime, end_date: datetime):
        self.start_date = pd.to_datetime(start_date, utc=True)
        self.end_date = pd.to_datetime(end_date, utc=True)

    @cache_dataframe()
    def so2(self, fuse: bool = False, smooth: bool = True) -> pd.DataFrame:
        """
        Retrieve SO2 gas flux for Whakaari. If fuse is True, the returned values
        will be a fusion of airborne and scandoas SO2 gas flux measurements. Otherwise,
        the returned values will be the airborne SO2 gas flux measurements from cospec
        if available and contouring otherwise.

        Parameters
        ----------
        fuse : bool
            If True, the returned values will be a fusion of airborne and scandoas
            SO2 gas flux measurements using a Kalman filter.
        smooth : bool
            If True, the returned values will be smoothed using a Kalman smoother.

        Returns
        -------
        pd.DataFrame
            Dataframe with SO2 gas flux values.
        """
        dataframes = {}

        tmp_df = []
        for method in ["cospec", "contouring", "flyspec", "mobile-doas"]:
            try:
                df = tilde_request(start_date=datetime(2000, 1, 1),
                                   end_date=self.end_date,
                                   domain="manualcollect",
                                   name="plume-SO2-gasflux",
                                   method=method,
                                   station="WI000", sensor="MC01")
            except ValueError:
                continue
            if len(df) > 0:
                tmp_df.append(df)
        if len(tmp_df) == 0:
            df = pd.DataFrame()
        else:
            df = tmp_df[0]
            if len(tmp_df) > 1:
                for df1 in tmp_df[1:]:
                    df = df.combine_first(df1)

        if not df.empty:
            dataframes['WI000'] = df

        if not fuse:
            return df * 86.4

        # Note: north-east point (WID01) is reporting consistently too low values
        for station in ["WID01", "WID02"]:
            try:
                df_ = tilde_request(start_date=datetime(2000, 1, 1),
                                    end_date=self.end_date,
                                    domain="scandoas",
                                    name="gasflux",
                                    method="reviewed",
                                    station=station,
                                    sensor="01")
            except ValueError:
                continue
            else:
                dataframes[station] = df_

        startdate = np.array([df.index.min()
                             for df in list(dataframes.values())]).min()
        enddate = np.array([df.index.max()
                           for df in list(dataframes.values())]).max()
        dates = pd.date_range(startdate, enddate, freq="1D")
        gasdata = {}
        interval = "1D"
        for key, df in dataframes.items():
            gasdata[key] = (
                df["obs"].groupby(pd.Grouper(freq=interval)
                                  ).median().reindex(dates).values
            )

        gasdata = pd.DataFrame(gasdata, index=dates)

        # define a covariance matrix for the observations that gives cospec values
        # a higher weight than mdoas values
        obs_cov = np.diag(np.array([5., 200., 20.])**2)
        model = SO2FusionModel(measurements=gasdata, initial_state=np.array([0]),
                               initial_cov=np.diag([[1]]),
                               obs_cov=obs_cov, k_states=1, k_posdef=1)
        # Set initial parameters for process noise
        initial_params = [0.5]
        # Fit the model
        result = model.fit(start_params=initial_params, disp=False)

        # Get the smoothed state estimates (filtered values)
        filtered_state = result.filtered_state[0]
        smoothed_state = result.smoothed_state[0]

        mean_trace = filtered_state
        error = np.sqrt(result.filtered_state_cov[0, 0])
        if smooth:
            mean_trace = smoothed_state
            error = np.sqrt(result.smoothed_state_cov[0, 0])
        so2df = pd.DataFrame(
            {"obs": mean_trace, "err": error}, index=gasdata.index)
        so2df.index.name = "dt"
        return so2df * 86.4

    @cache_dataframe()
    def co2(self) -> pd.DataFrame:
        df = tilde_request(start_date=datetime(2000, 1, 1),
                           end_date=self.end_date,
                           domain="manualcollect",
                           name="plume-CO2-gasflux",
                           method="contouring",
                           station="WI000", sensor="MC01")
        return df * 86.4

    @cache_dataframe()
    def h2s(self) -> pd.DataFrame:
        df = tilde_request(start_date=datetime(2000, 1, 1),
                           end_date=self.end_date,
                           domain="manualcollect",
                           name="plume-H2S-gasflux",
                           method="contouring",
                           station="WI000", sensor="MC01")
        return df * 86.4

class Seismicity(object):

    def __init__(self, start_date: datetime, end_date: datetime):
        self.start_date = pd.to_datetime(start_date, utc=True)
        self.end_date = pd.to_datetime(end_date, utc=True)

    def rsam(self) -> pd.DataFrame:
        """
        Load RSAM values NZ.WIZ.10.HHZ.The first time you call
        this it will download the data from Zenodo and store it
        locally.

        Returns:
        --------
            :param df: Dataframe with RSAM values.
            :type df: :class:`pandas.DataFrame`
        """
        rsam_fn = get_data('data/RSAM_NZ.WIZ.10.HHZ.csv')
        if not os.path.isfile(rsam_fn):
            # download the data from zenodo
            print('Downloading RSAM data from Zenodo')
            url = "https://zenodo.org/records/14759090/files/RSAM_NZ.WIZ.10.HHZ.csv"
            response = requests.get(url)
            if response.status_code == 200:
                with open(rsam_fn, "wb") as f:
                    f.write(response.content)
                print('Successfully downloaded RSAM data from Zenodo')
            else:
                print('Failed to download RSAM data from Zenodo')

        df = pd.read_csv(rsam_fn, skiprows=1, parse_dates=True, names=['dt', 'obs'], index_col=0,
                         date_format='ISO8601')
        if self.end_date.tzinfo is not None:
            df.index = df.index.tz_localize(timezone.utc)
        df = df[df.index <= str(self.end_date)]
        df = df[df.index >= str(self.start_date)]
        return df

    @cache_dataframe()
    def quakes(self):
        center_point = "177.186833+-37.523118"
        return wfs_request(self.start_date, self.end_date, radius=20000,
                           center_point=center_point, maxdepth=30)

class CraterLake(object):

    def __init__(self, start_date: datetime, end_date: datetime):
        self.start_date = pd.to_datetime(start_date, utc=True)
        self.end_date = pd.to_datetime(end_date, utc=True)

    @cache_dataframe()
    def temperature(self, resample: str = "D", interpolate: str = None,
                    exclude1995: bool = True, dropna: bool = True) -> pd.DataFrame:
        """
        Read crater lake temperature from Tilde (https://tilde.geonet.org.nz).

        Parameters:
        -----------
            :param resample: Average over resampling interval and
                            linearly interpolate in between. The
                            interval should be one of 'D', 'W',
                            or 'M'
            :type resample: str
            :param interpolate: Interpolate between resampled points.
                                Only takes effect if resample is not
                                None.
            :type interpolate: str

        Returns:
        --------
            :param df: A dataframe with the temperature time series.
            :type df: :class:`pandas.Dataframe`
        """
        dataframes = []
        # read in datalogger temperature data
        start_date = datetime(1954, 1, 1)
        sensors = [("envirosensor", "lake-temperature", "RU001", "-", "04"),
                   ("envirosensor", "lake-temperature", "RU001", "-", "01"),
                   ("envirosensor", "lake-temperature", "RU001", "snapshot", "06"),
                   ("envirosensor", "lake-temperature", "RU001", "snapshot", "07"),
                   ("manualcollect", "lake-temperature",
                    "RU001", "thermometer", "MC01"),
                   ("manualcollect", "lake-temperature", "RU001", "thermocouple", "MC01")]
        for domain, name, station, method, sensor in sensors:
            try:
                df = tilde_request(start_date=start_date, end_date=self.end_date,
                                   domain=domain, name=name,
                                   station=station, method=method, sensor=sensor)
            except ValueError:
                continue
            if len(df) > 0:
                dataframes.append(df)

        if len(dataframes) == 0:
            raise ValueError(
                f"No data found for lake temperature between {self.start_date} and {self.end_date}")

        df = dataframes[0]
        if len(dataframes) > 1:
            for df1 in dataframes[1:]:
                df = df.combine_first(df1)

        if resample is not None:
            df = df.resample("D").mean()
            if interpolate is not None:
                df = df.interpolate(interpolate)

        # remove observations following the 1995 eruption when
        # the measured temperatures do not represent a lake.
        if exclude1995:
            cond1 = df.index <= "1995-09-20 00:00:00"
            cond2 = df.index >= "2000-01-01 00:00:00"
            df = df[(cond1) | (cond2)]
        # get rid of nans
        if dropna:
            df.dropna(inplace=True)

        return df

    @cache_dataframe()
    def water_level(self, resample: str = "D", interpolate: str = None,
                    dropna: bool = True):
        """
        Read crater lake water level from Tilde (https://tilde.geonet.org.nz).
        """
        sensors = [("envirosensor", "lake-height", "RU001", "-", "02"),
                   ("envirosensor", "lake-height", "RU001", "-", "03")]
        dataframes = []
        for domain, name, station, method, sensor in sensors:
            try:
                df = tilde_request(start_date=datetime(2009, 4, 15, 2, 0, 0), end_date=self.end_date,
                                   domain=domain, name=name, station=station,
                                   method=method, sensor=sensor)
            except ValueError:
                continue
            if len(df) > 0:
                dataframes.append(df)
        try:
            df08 = tilde_request(start_date=datetime(2009, 4, 15, 2, 0, 0), end_date=self.end_date,
                                 domain="envirosensor", name="lake-height", station="RU001",
                                 method="snapshot", sensor="08")
        except ValueError:
            pass
        if len(df08) > 0:
            intercept, slope = lake_height_scaling()
            df08['obs'] = df08['obs'] * slope + intercept
            dataframes.append(df08)

        if len(dataframes) == 0:
            raise ValueError(
                f"No data found for lake height between {self.start_date} and {self.end_date}")
        df = dataframes[0]
        if len(dataframes) > 1:
            for df1 in dataframes[1:]:
                df = df.combine_first(df1)

        if resample is not None:
            df = df.resample("D").mean()
            if interpolate is not None:
                df = df.interpolate(interpolate)

        if dropna:
            df.dropna(inplace=True)

        return df

    @cache_dataframe()
    def water_analyte(self, analyte: str, resample: str = None, interpolate: str = None,
                      exclude1995: bool = True, drop_duplicates: bool = True) -> pd.DataFrame:
        """
        Download water analyte data from Tilde (https://tilde.geonet.org.nz).

        Parameters:
        -----------
            analyte : str
                The analyte to download. This can be one of
                    'Al, 'As', 'B', 'Br', 'Ca', 'Cl', 'Cs', 'F',
                    'Fe', 'H2S', 'K', 'Li', 'Mg', 'NH3', 'NO3-N', 'Na',
                    'PO4-P', 'Rb', 'SO4', 'SiO2', 'd18O', 'd2H', 'ph'
            resample : str
                Resample the data to a given interval. This can be anything
                allowed by :pandas.DataFrame.resample:. 
            interpolate : str
                Interpolate between points. This can be anything allowed by
                :pandas.DataFrame.interpolate:.
            exclude1995 : bool
                Exclude data from 1995 when there wasn't a lake.
        Returns:
        --------
            pandas.DataFrame
                A dataframe with the water analyte time series.
        """
        dataframes = []
        sensors = ["MC03", "MC04", "MC01"]
        for sensor in sensors:
            try:
                df = tilde_request(start_date=datetime(1964, 5, 9),
                                   end_date=self.end_date,
                                   domain="manualcollect",
                                   name=f"lake-{analyte}-conc",
                                   station="RU001",
                                   sensor=sensor)
            except ValueError:
                continue
            if len(df) > 0:
                dataframes.append(df)
        if len(dataframes) == 0:
            raise ValueError(
                f"No data found for {analyte} between {self.start_date} and {self.end_date}")
        df = dataframes[0]
        if len(dataframes) > 1:
            for df1 in dataframes[1:]:
                df = df.combine_first(df1)
        if resample is not None:
            df = df.resample(resample).mean()

        if interpolate is not None:
            df = df.interpolate(interpolate)

        if exclude1995:
            cond1 = df.index <= "1995-09-20 00:00:00"
            cond2 = df.index >= "2000-01-01 00:00:00"
            df = df[(cond1) | (cond2)]

        if drop_duplicates:
            # If there are multiple measurements on the same day, take the mean
            df = df.groupby(df.index.date).mean()
            df.index = pd.to_datetime(df.index, utc=True)
        return df

class Seismicity(object):

    def __init__(self, start_date: datetime, end_date: datetime):
        self.start_date = pd.to_datetime(start_date, utc=True)
        self.end_date = pd.to_datetime(end_date, utc=True)

    @cache_dataframe()
    def regional(self):
        """
        Load an earthquake catalogue for an area containing Waiouru and National Park.
        """

        polygon = (
            "175.32+-39.50,+175.32+-39.18,+175.77+-39.18,+175.77+-39.50,+175.32+-39.50"
        )
        tstart = datetime(1955, 1, 1, tzinfo=timezone.utc)
        return wfs_request(tstart, self.end_date, polygon=polygon)

    @cache_dataframe()
    def cone(self):
        """
        Load an earthquake catalogue for an area of 7 km around the summit.
        """
        radius = 7000
        center_point = "175.57+-39.28"
        tstart = datetime(1955, 1, 1, tzinfo=timezone.utc)
        return wfs_request(tstart, self.end_date, radius=radius, center_point=center_point)

    def rm_duplicates(self, cat1: pd.DataFrame, cat2: pd.DataFrame) -> pd.DataFrame:
        """
        Remove events from catalogue cat1 that are also in cat2 and return a new catalogue.

        Parameters:
        -----------
            cat1 : pandas.DataFrame
                The first catalogue.
            cat2 : pandas.DataFrame
                The second catalogue.

        Returns:
        --------
            pandas.DataFrame
                A new catalogue with the events from cat1 that are not in cat2.
        """
        idxs = []
        cat_tmp = cat1.copy()
        for i in range(cat_tmp.shape[0]):
            if cat_tmp.iloc[i].publicid in cat2.publicid.values:
                idxs.append(i)
        dates_to_drop = cat_tmp.index[idxs]
        cat_tmp.drop(dates_to_drop, inplace=True)
        return cat_tmp

    def daily_rsam(self, update=False) -> pd.DataFrame:
        """
        Load RSAM values from DRZ, MAVZ and FWVZ combined by scaling
        MAVZ and FWVZ RSAM values with DRZ.

        Returns:
        --------
            :param df: Dataframe with RSAM values.
            :type df: :class:`pandas.DataFrame`
        """
        df = pd.read_csv(get_data("data/ruapehu_rsam.csv"),
                         parse_dates=True, index_col=0)
        try:
            df.index = df.index.tz_localize(timezone.utc)
        except TypeError:
            pass
        if self.end_date > df.index[-1] and update:
            df_drz = pd.read_csv(
                get_data("data/DRZ_scaled_MAVZ.csv"), parse_dates=True, index_col=0
            )
            df_drz.rename(columns={"RSAM": "obs"}, inplace=True)
            url = "https://gessp.gns.cri.nz/vumt-api/feature?group=vumt&name=rsam"
            url += "&starttime=2007-01-01T00:00:00"
            url += f"&endtime={self.end_date.isoformat().split('+')[0]}"
            url += "&subdir=Ruapehu&subdir=FWVZ&subdir=HHZ"
            try:
                df_fwvz = pd.read_csv(
                    url,
                    parse_dates=True,
                    index_col=0,
                    date_format="ISO8601",
                )
            except Exception as e:
                print(url)
                raise e
            df_fwvz.rename(columns={"feature": "rsam"}, inplace=True)
            df_fwvz["rsam"] = np.where(
                df_fwvz["rsam"] > 1e30, np.nan, df_fwvz["rsam"])
            df_fwvz_daily = df_fwvz.resample("1D").mean()
            df_fwvz_daily_scaled = -2.3945 + 3.5062 * df_fwvz_daily
            df_fwvz_daily_scaled.rename(columns={"rsam": "obs"}, inplace=True)
            df_drz.index = df_drz.index.tz_localize(timezone.utc)
            df = df_drz.combine_first(df_fwvz_daily_scaled)
            # remove duplicated dates:
            df = df.loc[~df.index.duplicated(), :]
            df.to_csv(get_data("data/ruapehu_rsam.csv"))
        df = df[df.index <= str(self.end_date)]
        df = df[df.index >= str(self.start_date)]
        return df