Merge pull request #100 from ArtesiaWater/dev

Bug fix for BRO

hydropandas/extensions/geo.py

@@ -418,14 +418,11 @@ def get_lat_lon(self, in_epsg="epsg:28992", out_epsg="epsg:4326"):
         lon, lat : longitude and lattitude of x, y coordinates
         """
 
-        from pyproj import Proj, transform
+        from pyproj import Transformer
 
-        inProj = Proj(in_epsg)
-        outProj = Proj(out_epsg)
+        transformer = Transformer.from_crs(in_epsg, out_epsg)
 
         if np.isnan(self._obj.x) or np.isnan(self._obj.y):
-            lat, lon = np.nan, np.nan
-        else:
-            lat, lon = transform(inProj, outProj, self._obj.x, self._obj.y)
+            return np.nan, np.nan
 
-        return lat, lon
+        return transformer.transform(self._obj.x, self._obj.y)

hydropandas/extensions/gwobs.py

@@ -164,32 +164,32 @@ def get_modellayer_from_screen_depth(ftop, fbot, zvec, left=-999, right=999):
 
     # Piezometer in layer in which majority of screen is located
     if lay_fbot == lay_ftop:
-        logger.info(f"- selected layer {lay_fbot}")
+        logger.debug(f"- selected layer {lay_fbot}")
         return lay_fbot
 
     else:
         if lay_fbot == left and lay_ftop == right:
-            logger.info("- tube screen spans all layers. " "return nan")
+            logger.debug("- tube screen spans all layers. " "return nan")
             return np.nan
         elif lay_ftop == right:
-            logger.info(
+            logger.debug(
                 "- tube screen top higher than top layer. " f"selected layer {lay_fbot}"
             )
             return lay_fbot
 
         elif lay_fbot == left:
-            logger.info(
+            logger.debug(
                 "- tube screen bot lower than bottom layer. "
                 f"selected layer {lay_ftop}"
             )
             return lay_ftop
 
-        logger.info(
+        logger.debug(
             "- tube screen crosses layer boundary:\n"
             f"  - layers: {lay_ftop}, {lay_fbot}"
         )
 
-        logger.info(
+        logger.debug(
             f"- tube screen spans {lay_fbot - lay_ftop +1} layers."
             " checking length per layer\n"
             "  - length per layer:"
@@ -205,12 +205,12 @@ def get_modellayer_from_screen_depth(ftop, fbot, zvec, left=-999, right=999):
             else:
                 length_layers[i] = zvec[lay_ftop + i] - zvec[lay_ftop + 1 + i]
 
-            logger.info(f"    - lay {lay_ftop+i}: {length_layers[i]:.2f}")
+            logger.debug(f"    - lay {lay_ftop+i}: {length_layers[i]:.2f}")
 
         # choose layer with biggest length
         rel_layer = np.argmax(length_layers)
         lay_out = lay_ftop + rel_layer
-        logger.info(f"  - selected layer: {lay_out}")
+        logger.debug(f"  - selected layer: {lay_out}")
         return lay_out
 
     return np.nan
@@ -228,8 +228,8 @@ def get_zvec(x, y, gwf=None, ds=None):
     gwf : flopy.mf6.modflow.mfgwf.ModflowGwf
         modflow model with top and bottoms
     ds : xarray.Dataset
-        xarray Dataset with with top and bottoms, must have
-        dimensions 'x' and 'y' and variables 'top' and 'bot'
+        xarray Dataset typically create in nlmod. Must have
+        dimensions 'x' and 'y' and variables 'top' and 'botm'.
 
     Raises
     ------
@@ -277,10 +277,20 @@ def get_zvec(x, y, gwf=None, ds=None):
         xmin, xmax, ymin, ymax = ds.attrs["extent"]
         if (x < xmin) or (x > xmax) or (y < ymin) or (y > ymax):
             zvec = np.nan
+        elif ds.gridtype == "vertex":
+            import nlmod
+
+            cid = nlmod.dims.xy_to_icell2d((x, y), ds)
+            sel = ds.sel(icell2d=cid)
+            zvec = np.concatenate(([sel["top"].data], sel["botm"].data))
+            mask = np.isnan(zvec)
+            idx = np.where(~mask, np.arange(mask.size), 0)
+            np.maximum.accumulate(idx, axis=0, out=idx)
+            zvec[mask] = zvec[idx[mask]]
         else:
             sel = ds.sel(x=x, y=y, method="nearest")
             first_notna = np.nonzero(np.isfinite(sel["top"].data))[0][0]
-            zvec = np.concatenate([sel["top"].data[[first_notna]], sel["bot"].data])
+            zvec = np.concatenate([sel["top"].data[[first_notna]], sel["botm"].data])
             mask = np.isnan(zvec)
             idx = np.where(~mask, np.arange(mask.size), 0)
             np.maximum.accumulate(idx, axis=0, out=idx)
@@ -497,7 +507,7 @@ def get_modellayers(self, gwf=None, ds=None):
         """
         modellayers = []
         for o in self._obj.obs.values:
-            logger.info("-" * 10 + f"\n {o.name}:")
+            logger.debug("-" * 10 + f"\n {o.name}:")
 
             modellayers.append(o.gwobs.get_modellayer_modflow(gwf=gwf, ds=ds))
 
@@ -595,7 +605,11 @@ def get_regis_layer(self):
 
         # get index of regis model layer
         layer_i = get_modellayer_from_screen_depth(
-            self._obj.screen_top, self._obj.screen_bottom, zvec, left=-999, right=999,
+            self._obj.screen_top,
+            self._obj.screen_bottom,
+            zvec,
+            left=-999,
+            right=999,
         )
 
         if layer_i == 999:

hydropandas/extensions/plots.py

@@ -22,7 +22,12 @@ def __init__(self, oc_obj):
         self._obj = oc_obj
 
     def interactive_plots(
-        self, savedir, tmin=None, tmax=None, per_monitoring_well=True, **kwargs
+        self,
+        savedir="figures",
+        tmin=None,
+        tmax=None,
+        per_monitoring_well=True,
+        **kwargs,
     ):
         """Create interactive plots of the observations using bokeh.
 
@@ -41,7 +46,7 @@ def interactive_plots(
             will be passed to the Obs.to_interactive_plot method, options
             include:
 
-            - plot_columns : list of str
+            - cols : list of str or None
             - hoover_names : list of str
             - plot_freq : list of str
             - plot_legend_names : list of str
@@ -105,14 +110,14 @@ def interactive_plots(
                         return_filename=False,
                         **kwargs,
                     )
-                    logger.info(f"created iplot -> {o.name}")
+                    logger.debug(f"created iplot -> {o.name}")
                 except ValueError:
                     logger.error(f"{o.name} has no data between {tmin} and {tmax}")
                     o.iplot_fname = None
 
     def interactive_map(
         self,
-        plot_dir,
+        plot_dir="figures",
         m=None,
         tiles="OpenStreetMap",
         fname=None,
@@ -178,7 +183,7 @@ def interactive_map(
         **kwargs :
             will be passed to the to_interactive_plots method options are:
 
-            - plot_columns : list of str
+            - cols : list of str or None
             - hoover_names : list of str
             - plot_legend_names : list of str
             - plot_freq : list of str

hydropandas/extensions/stats.py

@@ -322,8 +322,8 @@ def get_seasonal_stat(
         return df
 
     def obs_per_year(self, col=None):
-        """get number of observations per year 
-        
+        """get number of observations per year
+
         Parameters
         ----------
         col : str or None, optional
@@ -361,7 +361,8 @@ def consecutive_obs_years(obs_per_year, min_obs=12):
         obs_per_year_all = pd.Series(index=[dt.datetime.now().year], data=0)
     else:
         obs_per_year_all = pd.Series(
-            dtype=float, index=range(obs_per_year.index[0], obs_per_year.index[-1] + 1),
+            dtype=float,
+            index=range(obs_per_year.index[0], obs_per_year.index[-1] + 1),
         )
         obs_per_year_all.loc[obs_per_year.index] = obs_per_year
 

hydropandas/io/io_bro.py

@@ -91,9 +91,9 @@ class of the observations, so far only GroundwaterObs is supported
 
 
 def get_bro_groundwater(bro_id, tube_nr=None, only_metadata=False, **kwargs):
-    """ get bro groundwater measurement from a GLD id or a GMW id with a
-    filter number. 
-    
+    """get bro groundwater measurement from a GLD id or a GMW id with a
+    filter number.
+
 
     Parameters
     ----------
@@ -149,7 +149,7 @@ def get_bro_groundwater(bro_id, tube_nr=None, only_metadata=False, **kwargs):
 
 
 def get_gld_id_from_gmw(bro_id, tube_nr):
-    """ get bro_id of a grondwterstandendossier (gld) from a bro_id of a 
+    """get bro_id of a grondwterstandendossier (gld) from a bro_id of a
     grondwatermonitoringsput (gmw).
 
     Parameters
@@ -205,9 +205,9 @@ def get_gld_id_from_gmw(bro_id, tube_nr):
 def measurements_from_gld(
     bro_id, tmin=None, tmax=None, to_wintertime=True, drop_duplicate_times=True
 ):
-    """ get measurements and metadata from a grondwaterstandonderzoek (gld)
+    """get measurements and metadata from a grondwaterstandonderzoek (gld)
     bro_id
-    
+
 
     Parameters
     ----------
@@ -265,6 +265,7 @@ def measurements_from_gld(
     }
 
     tree = xml.etree.ElementTree.fromstring(req.text)
+
     glds = tree.findall(".//ns11:GLD_O", ns)
     if len(glds) != 1:
         raise (Exception("Only one gld supported"))
@@ -275,19 +276,25 @@ def measurements_from_gld(
     meta["tube_nr"] = int(
         gld.find("ns11:monitoringPoint//gldcommon:tubeNumber", ns).text
     )
-    meta["monitoringsnet"] = gld.find(
-        "ns11:groundwaterMonitoringNet//gldcommon:broId", ns
-    ).text
+    gmn = gld.find("ns11:groundwaterMonitoringNet//gldcommon:broId", ns)
+    if gmn is None:
+        meta["monitoringsnet"] = None
+    else:
+        meta["monitoringsnet"] = gmn.text
 
     # get observations
     msts = "ns11:observation//om:result//waterml:MeasurementTimeseries"
     times = [time.text for time in gld.findall(f"{msts}//waterml:time", ns)]
-    values = [float(value.text) for value in gld.findall(f"{msts}//waterml:value", ns)]
-    qualifiers = [q.text for q in gld.findall(f"{msts}//swe:value", ns)]
+    values = [
+        np.nan if value.text is None else float(value.text)
+        for value in gld.findall(f"{msts}//waterml:value", ns)
+    ]
+    qualifiers = [q.text for q in gld.findall(f"{msts}//swe:Category//swe:value", ns)]
 
     # to dataframe
     df = pd.DataFrame(
-        index=pd.to_datetime(times), data={"values": values, "qualifier": qualifiers},
+        index=pd.to_datetime(times),
+        data={"values": values, "qualifier": qualifiers},
     )
 
     # wintertime
@@ -301,7 +308,7 @@ def measurements_from_gld(
     if df.index.has_duplicates and drop_duplicate_times:
         duplicates = df.index.duplicated(keep="first")
         message = "{} contains {} duplicates (of {}). Keeping only first values."
-        logger.warning(message.format(bro_id, duplicates.sum(), len(df)))
+        logger.info(message.format(bro_id, duplicates.sum(), len(df)))
         df = df[~duplicates]
 
     df = df.sort_index()
@@ -316,8 +323,8 @@ def measurements_from_gld(
 
 
 def get_full_metadata_from_gmw(bro_id, tube_nr):
-    """ get metadata for a groundwater monitoring well.
-    
+    """get metadata for a groundwater monitoring well.
+
 
     Parameters
     ----------
@@ -406,9 +413,9 @@ def get_full_metadata_from_gmw(bro_id, tube_nr):
 
 
 def get_metadata_from_gmw(bro_id, tube_nr):
-    """ get selection of metadata for a groundwater monitoring well.
+    """get selection of metadata for a groundwater monitoring well.
     coordinates, ground_level, tube_top and tube screen
-    
+
 
     Parameters
     ----------
@@ -508,8 +515,8 @@ def get_obs_list_from_extent(
     epsg=28992,
     ignore_max_obs=False,
 ):
-    """ get a list of gmw observations within an extent.
-    
+    """get a list of gmw observations within an extent.
+
 
     Parameters
     ----------
@@ -534,7 +541,7 @@ class of the observations, e.g. GroundwaterObs or WaterlvlObs
 
     Raises
     ------
-    
+
         DESCRIPTION.
 
     Returns
@@ -580,22 +587,24 @@ class of the observations, e.g. GroundwaterObs or WaterlvlObs
     if tree.find(".//brocom:responseType", ns).text == "rejection":
         raise RuntimeError(tree.find(".//brocom:rejectionReason", ns).text)
 
-    gmws = tree.findall(".//dsgmw:GMW_C", ns)
+    gmws_ids = np.unique(
+        [gmw.text for gmw in tree.findall(".//dsgmw:GMW_C//brocom:broId", ns)]
+    )
 
-    if len(gmws) > 1000:
+    if len(gmws_ids) > 1000:
         ans = input(
-            f"You requested to download {len(gmws)} observations, this can take a while. Are you sure you want to continue [Y/n]? "
+            f"You requested to download {len(gmws_ids)} observations, this can take a while. Are you sure you want to continue [Y/n]? "
         )
         if ans not in ["Y", "y", "yes", "Yes", "YES"]:
             return []
 
     obs_list = []
-    gmw_list = []
-    for gmw in tqdm(gmws):
-        gmw_id = gmw.find("brocom:broId", ns).text
-        if gmw_id in gmw_list:
-            logger.debug(f"gmw_id {gmw_id} already read, skipping")
-            continue
+    for gmw_id in tqdm(gmws_ids):
+        gmws = tree.findall(f'.//*[brocom:broId="{gmw_id}"]', ns)
+        if len(gmws) < 1:
+            raise RuntimeError("unexpected")
+        else:
+            gmw = gmws[0]
 
         ntubes = int(gmw.find("dsgmw:numberOfMonitoringTubes", ns).text)
         for tube_nr in range(1, ntubes + 1):
@@ -614,6 +623,5 @@ class of the observations, e.g. GroundwaterObs or WaterlvlObs
                     obs_list.append(o)
             else:
                 obs_list.append(o)
-        gmw_list.append(gmw_id)
 
     return obs_list