Merge pull request #65 from ArtesiaWater/dev

minor fixes for knmi meteo var and plots

hydropandas/extensions/plots.py

@@ -43,10 +43,12 @@ def interactive_plots(self, savedir,
 
             - plot_columns : list of str
             - hoover_names : list of str
-            - plot_freq : str
-            - plot_legend_name : str
+            - plot_freq : list of str
+            - plot_legend_names : list of str
+            - markers : list of str
+            - hoover_names : list of str
+            - plot_colors : list of str
             - ylabel : str
-            - color : str
             - add_filter_to_legend : boolean
         """
         _color_cycle = (
@@ -84,7 +86,7 @@ def interactive_plots(self, savedir,
                     p = o.plots.interactive_plot(savedir=savedir,
                                                  p=p,
                                                  tmin=tmin, tmax=tmax,
-                                                 colors=[_color_cycle[i + 1]],
+                                                 plot_colors=[_color_cycle[i + 1]],
                                                  return_filename=False,
                                                  **kwargs)
                     logger.info(f'created iplot -> {o.name}')
@@ -157,10 +159,13 @@ def interactive_map(self, plot_dir, m=None,
 
             - plot_columns : list of str
             - hoover_names : list of str
-            - plot_legend_name : str
+            - plot_legend_names : list of str
+            - plot_freq : list of str
+            - markers : list of str
+            - hoover_names : list of str
+            - plot_colors : list of str
             - ylabel : str
             - add_filter_to_legend : boolean
-            - plot_freq : str
             - tmin : dt.datetime
             - tmax : dt.datetime
 
@@ -272,7 +277,7 @@ def interactive_plot(self,
                          hoover_names=['Peil'],
                          hoover_date_format="%Y-%m-%d",
                          ylabel='m NAP',
-                         colors=['blue'],
+                         plot_colors=['blue'],
                          add_filter_to_legend=False,
                          return_filename=False):
         """Create an interactive plot of the observations using bokeh.
@@ -307,8 +312,8 @@ def interactive_plot(self,
             date format to use when hoovering over a plot
         ylabel : str, optional
             label on the y-axis
-        colors : list of str, optional
-            colors used for the plots
+        plot_colors : list of str, optional
+            plot_colors used for the plots
         add_filter_to_legend : boolean, optional
             if True the attributes bovenkant_filter and onderkant_filter
             are added to the legend name
@@ -345,8 +350,13 @@ def interactive_plot(self,
             xcol = 'index'
 
         # get color
-        if len(colors) < len(plot_columns):
-            colors = colors * len(plot_columns)
+        if len(plot_colors) < len(plot_columns):
+            plot_colors = plot_colors * len(plot_columns)
+
+        # get base for hoover tooltips
+        plots = []
+        tooltips = []
+        tooltips.append(('date', "@date"))
 
         # plot multiple columns
         for i, column in enumerate(plot_columns):
@@ -366,26 +376,24 @@ def interactive_plot(self,
                     plot_df[[column, 'date']].resample(plot_freq[i]).first())
 
             # plot data
+
             if markers[i] == 'line':
-                p.line(xcol, column, source=source, color=colors[i],
-                       legend_label=lname,
-                       alpha=0.8, muted_alpha=0.2)
+                plots.append(p.line(xcol, column, source=source, color=plot_colors[i],
+                               legend_label=lname,
+                               alpha=0.8, muted_alpha=0.2))
             elif markers[i] == 'circle':
-                p.circle(xcol, column, source=source, color=colors[i],
-                         legend_label=lname,
-                         alpha=0.8, muted_alpha=0.2)
+                plots.append(p.circle(xcol, column, source=source, color=plot_colors[i],
+                                 legend_label=lname,
+                                 alpha=0.8, muted_alpha=0.2))
             else:
                 raise NotImplementedError("marker '{}' invalid. Only line and"
                                           "circle are currently available".format(markers[i]))
 
-        # hoover options
-        tooltips = []
-        for i, column in enumerate(plot_columns):
-            tooltips.append((hoover_names[i], "@{}".format(column)))
-        tooltips.append(('date', "@date"))
-        hover = HoverTool(tooltips=tooltips)
-
-        p.add_tools(hover)
+            # add columns to hoover tooltips
+            tooltips_p = tooltips.copy()
+            tooltips_p.append((hoover_names[i], "@{}".format(column)))
+            hover = HoverTool(renderers=[plots[i]], tooltips=tooltips_p, mode='vline')
+            p.add_tools(hover)
 
         p.legend.location = "top_left"
         p.legend.click_policy = "mute"

hydropandas/io/io_knmi.py

@@ -275,13 +275,15 @@ def _check_latest_measurement_date_RD_debilt(use_api=True):
     knmi_df = knmi_df.dropna()
     if knmi_df.empty:
         raise ValueError(
-            f'knmi station de Bilt has no RD measurements in the past {look_back_days} days.')
+            'knmi station de Bilt has no RD measurements '
+            f'in the past {look_back_days} days.')
 
     last_measurement_date_debilt = knmi_df.index[-1]
 
     logger.info(f'last RD measurement available at the Bilt is from'
                 f' {last_measurement_date_debilt.strftime("%Y-%m-%d")}')
-    logger.info('assuming no measurements are available at other stations before this date')
+    logger.info('assuming no measurements are available at '
+                'other stations before this date')
 
     return last_measurement_date_debilt
 
@@ -349,7 +351,9 @@ def download_knmi_data(stn, stn_name=None,
     # convert possible integer to string
     stn = str(stn)
 
-    logger.info(f'download knmi {meteo_var} data from station {stn}-{stn_name} between {start} and {end}')
+    logger.info(
+        f'download knmi {meteo_var} data from station '
+        f'{stn}-{stn_name} between {start} and {end}')
 
     # define variables
     knmi_df = pd.DataFrame()
@@ -361,7 +365,8 @@ def download_knmi_data(stn, stn_name=None,
         if use_api:
             if interval.startswith('hour'):
                 # hourly data from meteorological stations
-                knmi_df = get_knmi_hourly_api(URL_HOURLY_METEO, stn, meteo_var, start, end)
+                knmi_df = get_knmi_hourly_api(
+                    URL_HOURLY_METEO, stn, meteo_var, start, end)
 
             elif meteo_var == 'RD':
                 # daily data from rainfall-stations
@@ -949,27 +954,11 @@ def get_knmi_timeseries_xy(x, y, meteo_var, start, end,
     stations = get_stations(meteo_var=meteo_var)
     stn = get_nearest_stations_xy(x, y, meteo_var, stations=stations)[0]
 
-    # download data
-    if fill_missing_obs:
-        knmi_df, variables, station_meta = \
-            fill_missing_measurements(stn, stn_name, meteo_var, start, end,
-                                      interval, raise_exceptions)
-    else:
-        knmi_df, variables, station_meta = \
-            download_knmi_data(stn, stn_name, meteo_var, start, end,
-                               interval, inseason, raise_exceptions)
-
-    if not station_meta is None:
-        meta = station_meta.to_dict()
-    else:
-        meta = {}
-    meta.update(variables)
-
-    # set metadata
-    name = meteo_var + ' ' + stations.loc[stn, 'naam']
-    x = stations.loc[stn, 'x']
-    y = stations.loc[stn, 'y']
-    meta.update({'x': x, 'y': y, 'station': stn, 'name': name})
+    knmi_df, meta = get_knmi_timeseries_stn(stn, meteo_var, start, end,
+                                            fill_missing_obs=fill_missing_obs,
+                                            interval=interval,
+                                            inseason=inseason,
+                                            raise_exceptions=raise_exceptions)
 
     return knmi_df, meta
 
@@ -1036,7 +1025,12 @@ def get_knmi_timeseries_stn(stn, meteo_var, start, end,
     # set metadata
     x = stations.loc[stn, 'x']
     y = stations.loc[stn, 'y']
-    meta.update({'x': x, 'y': y, 'station': stn, 'name': stn_name})
+    meta.update({'x': x,
+                 'y': y,
+                 'station': stn,
+                 'name': f"{meteo_var}_{stn_name}",
+                 "variable": meteo_var
+                 })
 
     return knmi_df, meta
 
@@ -1176,7 +1170,8 @@ def add_missing_indices(knmi_df, stn, start, end):
                                  day=start.day, hour=knmi_df.index[0].hour,
                                  minute=knmi_df.index[0].minute,
                                  second=knmi_df.index[0].second)
-        logger.info(f'station {stn} has no measurements before {knmi_df.index[0]}')
+        logger.info(
+            f'station {stn} has no measurements before {knmi_df.index[0]}')
 
     if (end - knmi_df.index[-1]).days < 2:
         new_end = knmi_df.index[-1]
@@ -1185,7 +1180,8 @@ def add_missing_indices(knmi_df, stn, start, end):
                                hour=knmi_df.index[-1].hour,
                                minute=knmi_df.index[-1].minute,
                                second=knmi_df.index[-1].second)
-        logger.info(f'station {stn} has no measurements after {knmi_df.index[-1]}')
+        logger.info(
+            f'station {stn} has no measurements after {knmi_df.index[-1]}')
 
     # add missing indices
     new_index = pd.date_range(new_start, new_end, freq='D')
@@ -1255,7 +1251,8 @@ def fill_missing_measurements(stn, stn_name=None, meteo_var='RD',
     # if the first station cannot be read, read another station as the first
     ignore = [stn]
     while knmi_df.empty:
-        logger.info(f'station {stn} has no measurements between {start} and {end}')
+        logger.info(
+            f'station {stn} has no measurements between {start} and {end}')
         logger.info('trying to get measurements from nearest station')
         stn = get_nearest_station_df(
             stations.loc[[stn]], meteo_var=meteo_var, ignore=ignore)[0]
@@ -1280,11 +1277,11 @@ def fill_missing_measurements(stn, stn_name=None, meteo_var='RD',
             stations.loc[[stn]], meteo_var=meteo_var, ignore=ignore)
 
         logger.info(f'trying to fill {missing.sum()} '
-                     f'measurements with station {stn_comp}')
+                    f'measurements with station {stn_comp}')
 
         if stn_comp is None:
             logger.info('could not fill all missing measurements there are '
-                         'no stations left to check')
+                        'no stations left to check')
 
             missing[:] = False
             break

hydropandas/io/io_modflow.py

@@ -2,6 +2,7 @@
 import warnings
 
 import numpy as np
+from scipy.interpolate import griddata
 
 from .. import util
 from ..observation import ModelObs
@@ -84,7 +85,8 @@ def read_imod_results(obs_collection, ml, runfile, mtime, model_ws,
 
 
 def read_modflow_results(obs_collection, ml, hds_arr, mtime,
-                         modelname='', nlay=None, exclude_layers=None):
+                         modelname='', nlay=None, exclude_layers=None,
+                         method="linear"):
     """Read modflow groundwater heads at points in obs_collection.
 
     Parameters
@@ -103,7 +105,9 @@ def read_modflow_results(obs_collection, ml, hds_arr, mtime,
         number of layers if None the number of layers from ml is used.
     exclude_layers : list of int, optional
         exclude the observations in these modellayers
-
+    method : str, optional
+        interpolation method, either 'linear' or 'nearest',
+        default is linear.
     """
     if ml.modelgrid.grid_type == 'structured':
         if ml.modelgrid.xoffset == 0 or ml.modelgrid.yoffset == 0:
@@ -124,7 +128,8 @@ def read_modflow_results(obs_collection, ml, hds_arr, mtime,
         xy = np.array([xmid, ymid]).T
 
     uv = obs_collection.loc[:, ("x", "y")].dropna(how="any", axis=0).values
-    vtx, wts = util.interp_weights(xy, uv)
+    if method == "linear":
+        vtx, wts = util.interp_weights(xy, uv)
 
     # get interpolated timeseries from hds_arr
     hm_ts = np.nan * np.ones((obs_collection.shape[0], hds_arr.shape[0]))
@@ -141,8 +146,14 @@ def read_modflow_results(obs_collection, ml, hds_arr, mtime,
         for t in range(hds_arr.shape[0]):
             ihds = hds_arr[t, m]
             ihds[ihds <= -999.] = np.nan
-            hm = util.interpolate(ihds, vtx, wts)
-            hm_ts[mask, t] = hm[mask]
+            if method == "linear":
+                hm = util.interpolate(ihds, vtx, wts)
+                hm_ts[mask, t] = hm[mask]
+            elif method == "nearest":
+                hm_ts[mask, t] = griddata(
+                    xy, ihds.ravel(), uv, method=method)[mask]
+            else:
+                raise ValueError(f"Unknown method: '{method}'")
 
     mo_list = []
     for i, name in enumerate(obs_collection.index):

hydropandas/observation.py

@@ -585,11 +585,12 @@ class KnmiObs(Obs):
     Subclass of the Obs class
     """
 
-    _metadata = Obs._metadata + ['station']
+    _metadata = Obs._metadata + ['station', 'meteo_var']
 
     def __init__(self, *args, **kwargs):
 
         self.station = kwargs.pop('station', np.nan)
+        self.meteo_var = kwargs.pop('meteo_var', '')
 
         super(KnmiObs, self).__init__(*args, **kwargs)
 
@@ -598,7 +599,7 @@ def _constructor(self):
         return KnmiObs
 
     @classmethod
-    def from_knmi(cls, stn, variable, startdate=None, enddate=None,
+    def from_knmi(cls, stn, meteo_var, startdate=None, enddate=None,
                   fill_missing_obs=True, interval='daily', inseason=False,
                   use_api=True, raise_exceptions=True):
         """Get a KnmiObs object.
@@ -607,7 +608,7 @@ def from_knmi(cls, stn, variable, startdate=None, enddate=None,
         ----------
         stn : int or str
             measurement station e.g. 829.
-        variable : str, optional
+        meteo_var : str, optional
             observation type e.g. "RD" or "EV24". See list with all variables
             below.
         startdate : str, datetime or None, optional
@@ -677,16 +678,16 @@ def from_knmi(cls, stn, variable, startdate=None, enddate=None,
         from .io import io_knmi
 
         ts, meta = io_knmi.get_knmi_timeseries_stn(
-            stn, variable, startdate, enddate, fill_missing_obs,
+            stn, meteo_var, startdate, enddate, fill_missing_obs,
             interval=interval, inseason=inseason,
             use_api=use_api,
             raise_exceptions=raise_exceptions)
 
         return cls(ts, meta=meta, station=meta['station'], x=meta['x'],
-                   y=meta['y'], name=meta['name'])
+                   y=meta['y'], name=meta['name'], meteo_var=meteo_var)
 
     @classmethod
-    def from_nearest_xy(cls, x, y, variable, startdate=None, enddate=None,
+    def from_nearest_xy(cls, x, y, meteo_var, startdate=None, enddate=None,
                         fill_missing_obs=True, interval='daily',
                         inseason=False, raise_exceptions=False):
         """Get KnmiObs object with measurements from station closest to the
@@ -698,7 +699,7 @@ def from_nearest_xy(cls, x, y, variable, startdate=None, enddate=None,
             x coördinate in m RD.
         y : int or float
             y coördinate in m RD.
-        variable : str
+        meteo_var : str
             e.g. 'EV24'.
         startdate : str, datetime or None, optional
             start date of observations. The default is None.
@@ -721,15 +722,16 @@ def from_nearest_xy(cls, x, y, variable, startdate=None, enddate=None,
         from .io import io_knmi
 
         ts, meta = io_knmi.get_knmi_timeseries_xy(
-            x, y, variable, startdate, enddate, fill_missing_obs,
+            x, y, meteo_var, startdate, enddate, stn_name=None,
+            fill_missing_obs=fill_missing_obs,
             interval=interval, inseason=inseason,
             raise_exceptions=raise_exceptions)
 
         return cls(ts, meta=meta, station=meta['station'], x=meta['x'],
-                   y=meta['y'], name=meta['name'])
+                   y=meta['y'], name=meta['name'], meteo_var=meteo_var)
 
     @classmethod
-    def from_obs(cls, obs, variable, startdate=None, enddate=None,
+    def from_obs(cls, obs, meteo_var, startdate=None, enddate=None,
                  fill_missing_obs=True, interval='daily', inseason=False,
                  raise_exceptions=False):
 
@@ -744,9 +746,9 @@ def from_obs(cls, obs, variable, startdate=None, enddate=None,
             enddate = obs.index[-1]
 
         ts, meta = io_knmi.get_knmi_timeseries_xy(
-            x, y, variable, startdate, enddate, fill_missing_obs,
+            x, y, meteo_var, startdate, enddate, fill_missing_obs,
             interval=interval, inseason=inseason,
             raise_exceptions=raise_exceptions)
 
         return cls(ts, meta=meta, station=meta['station'], x=meta['x'],
-                   y=meta['y'], name=meta['name'])
+                   y=meta['y'], name=meta['name'], meteo_var=meteo_var)

hydropandas/version.py

@@ -1 +1 @@
-__version__ = '0.4.0'
+__version__ = '0.4.1'