Merge pull request #197 from ArtesiaWater/dev

Update to 0.6.1

docs/examples/06_gridding_vector_data.ipynb

@@ -1,6 +1,7 @@
 {
  "cells": [
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -19,27 +20,20 @@
    "outputs": [],
    "source": [
     "import nlmod\n",
-    "from nlmod import resample\n",
     "import numpy as np\n",
     "import xarray as xr\n",
-    "import flopy\n",
-    "import warnings\n",
     "\n",
-    "\n",
-    "from matplotlib.colors import Normalize\n",
-    "from matplotlib.patches import Polygon\n",
-    "from matplotlib.collections import PatchCollection\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
     "import geopandas as gpd\n",
     "from shapely.geometry import LineString, Point\n",
     "from shapely.geometry import Polygon as shp_polygon\n",
-    "from scipy.interpolate import RectBivariateSpline\n",
     "\n",
     "from IPython.display import display"
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -63,6 +57,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -121,13 +116,15 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "### Points"
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -168,6 +165,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -185,12 +183,12 @@
     "sim = nlmod.sim.sim(ds)\n",
     "gwf = nlmod.gwf.gwf(ds, sim)\n",
     "dis = nlmod.gwf.dis(ds, gwf)\n",
-    "da1 = nlmod.grid.gdf_to_data_array_struc(\n",
-    "    point_gdf, gwf, field=\"values\", interp_method=\"nearest\"\n",
-    ")\n",
-    "da2 = nlmod.grid.gdf_to_data_array_struc(\n",
-    "    point_gdf, gwf, field=\"values\", interp_method=\"linear\"\n",
+    "da1 = nlmod.grid.gdf_to_da(\n",
+    "    point_gdf, ds, column=\"values\", agg_method=\"nearest\"\n",
     ")\n",
+    "da2 = xr.DataArray(np.nan, dims=(\"y\", \"x\"), coords={\"y\": ds.y, \"x\": ds.x})\n",
+    "da2.values = nlmod.grid.interpolate_gdf_to_array(point_gdf, gwf, field='values', \n",
+    "method='linear')\n",
     "\n",
     "vmin = min(da1.min(), da2.min())\n",
     "vmax = max(da1.max(), da2.max())\n",
@@ -211,6 +209,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -251,6 +250,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -267,7 +267,7 @@
     "\n",
     "da1 = nlmod.grid.gdf_to_da(pol_gdf, ds, \"values\", agg_method=\"max_area\")\n",
     "da2 = nlmod.grid.gdf_to_da(pol_gdf, ds, \"values\", agg_method=\"area_weighted\")\n",
-    "da3 = nlmod.grid.gdf_to_data_array_struc(pol_gdf, gwf, \"values\", agg_method=\"nearest\")\n",
+    "da3 = nlmod.grid.gdf_to_da(pol_gdf, ds, \"values\", agg_method=\"nearest\")\n",
     "\n",
     "vmin = min(da1.min(), da2.min(), da3.min())\n",
     "vmax = max(da1.max(), da2.max(), da3.max())\n",
@@ -291,6 +291,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -328,6 +329,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -373,6 +375,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -416,6 +419,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -438,6 +442,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -461,6 +466,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -496,6 +502,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -517,6 +524,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -540,6 +548,7 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
@@ -566,8 +575,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "nlmod",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.4"
   }
  },
  "nbformat": 4,

nlmod/dims/base.py

@@ -319,11 +319,15 @@ def _get_structured_grid_ds(
         "y": ycenters,
         "layer": range(nlay),
     }
+
     if angrot != 0.0:
         affine = resample.get_affine_mod_to_world(attrs)
         xc, yc = affine * np.meshgrid(xcenters, ycenters)
         coords["xc"] = (("y", "x"), xc)
         coords["yc"] = (("y", "x"), yc)
+    else:
+        coords["x"] += xorigin
+        coords["y"] += yorigin
 
     dims = ("layer", "y", "x")
     ds = xr.Dataset(
@@ -439,7 +443,7 @@ def _get_vertex_grid_ds(
     else:
         layers = nlay
 
-    coords = {"x": x, "y": y, "layer": layers}
+    coords = {"layer": layers, "y": y, "x": x}
     dims = ("layer", "icell2d")
     ds = xr.Dataset(
         data_vars=dict(

nlmod/dims/grid.py

@@ -1058,13 +1058,22 @@ def gdf_to_da(
             raise ValueError(
                 "multiple geometries in one cell please define aggregation method"
             )
-        gdf_agg = aggregate_vector_per_cell(gdf_cellid, {column: agg_method})
+        if agg_method in ["nearest"]:
+            modelgrid = modelgrid_from_ds(ds)
+            gdf_agg = aggregate_vector_per_cell(
+                gdf_cellid, {column: agg_method}, modelgrid
+            )
+        else:
+            gdf_agg = aggregate_vector_per_cell(gdf_cellid, {column: agg_method})
     else:
         # aggregation not neccesary
         gdf_agg = gdf_cellid[[column]]
-        gdf_agg.set_index(
-            pd.MultiIndex.from_tuples(gdf_cellid.cellid.values), inplace=True
-        )
+        if isinstance(gdf_cellid.cellid.iloc[0], tuple):
+            gdf_agg.set_index(
+                pd.MultiIndex.from_tuples(gdf_cellid.cellid.values), inplace=True
+            )
+        else:
+            gdf_agg.set_index(gdf_cellid.cellid.values, inplace=True)
     da = util.get_da_from_da_ds(ds, dims=ds.top.dims, data=fill_value)
     for ind, row in gdf_agg.iterrows():
         da.values[ind] = row[column]
@@ -1135,17 +1144,22 @@ def _agg_length_weighted(gdf, col):
     return aw
 
 
-def _agg_nearest(gdf, col, gwf):
-    cid = gdf["cellid"].values[0]
-    cellcenter = Point(
-        gwf.modelgrid.xcellcenters[0][cid[1]],
-        gwf.modelgrid.ycellcenters[:, 0][cid[0]],
-    )
-    val = gdf.iloc[gdf.distance(cellcenter).argmin()].loc[col]
+def _agg_nearest(gdf, col, modelgrid):
+    if modelgrid.grid_type == "structured":
+        cid = gdf["cellid"].values[0]
+        cellcenter = Point(
+            modelgrid.xcellcenters[0, cid[1]], modelgrid.ycellcenters[cid[0], 0]
+        )
+        val = gdf.iloc[gdf.distance(cellcenter).argmin()].loc[col]
+    elif modelgrid.grid_type == "vertex":
+        cid = gdf["cellid"].values[0]
+        cellcenter = Point(modelgrid.xcellcenters[cid], modelgrid.ycellcenters[cid])
+        val = gdf.iloc[gdf.distance(cellcenter).argmin()].loc[col]
+
     return val
 
 
-def _get_aggregates_values(group, fields_methods, gwf=None):
+def _get_aggregates_values(group, fields_methods, modelgrid=None):
     agg_dic = {}
     for field, method in fields_methods.items():
         # aggregation is only necesary if group shape is greater than 1
@@ -1160,7 +1174,7 @@ def _get_aggregates_values(group, fields_methods, gwf=None):
         elif method == "sum":
             agg_dic[field] = group[field].sum()
         elif method == "nearest":
-            agg_dic[field] = _agg_nearest(group, field, gwf)
+            agg_dic[field] = _agg_nearest(group, field, modelgrid)
         elif method == "length_weighted":  # only for lines
             agg_dic[field] = _agg_length_weighted(group, field)
         elif method == "max_length":  # only for lines
@@ -1177,7 +1191,7 @@ def _get_aggregates_values(group, fields_methods, gwf=None):
     return agg_dic
 
 
-def aggregate_vector_per_cell(gdf, fields_methods, gwf=None):
+def aggregate_vector_per_cell(gdf, fields_methods, modelgrid=None):
     """Aggregate vector features per cell.
 
     Parameters
@@ -1189,7 +1203,7 @@ def aggregate_vector_per_cell(gdf, fields_methods, gwf=None):
         aggregation methods can be:
         max, min, mean, sum, length_weighted (lines), max_length (lines),
         area_weighted (polygon), max_area (polygon).
-    gwf : flopy Groundwater flow model
+    modelgrid : flopy Groundwater flow modelgrid
         only necesary if one of the field methods is 'nearest'
 
     Returns
@@ -1227,7 +1241,7 @@ def aggregate_vector_per_cell(gdf, fields_methods, gwf=None):
     gr = gdf.groupby(by="cellid")
     celldata = pd.DataFrame(index=gr.groups.keys())
     for cid, group in tqdm(gr, desc="Aggregate vector data"):
-        agg_dic = _get_aggregates_values(group, fields_methods, gwf)
+        agg_dic = _get_aggregates_values(group, fields_methods, modelgrid)
         for key, item in agg_dic.items():
             celldata.loc[cid, key] = item
 

nlmod/version.py

@@ -1,7 +1,7 @@
 from importlib import metadata
 from platform import python_version
 
-__version__ = "0.6.0"
+__version__ = "0.6.1"
 
 
 def show_versions() -> None:

tests/test_008_waterschappen.py

@@ -5,7 +5,7 @@
 """
 
 import pytest
-
+import matplotlib
 import nlmod
 
 
@@ -48,13 +48,23 @@ def test_download_peilgebieden(plot=True):
 
     if plot:
         # plot the winter_stage
-        ax = waterboards.plot(edgecolor="k", facecolor="none")
+        f, ax = nlmod.plot.get_map([9000, 279000, 304000, 623000], base=100000)
+        waterboards.plot(edgecolor="k", facecolor="none", ax=ax)
+        norm = matplotlib.colors.Normalize(-10.0, 20.0)
+        cmap = "viridis"
         for wb in waterboards.index:
             if wb in gdf:
                 # gdf[wb].plot(ax=ax, zorder=0)
-                gdf[wb].plot("winter_stage", ax=ax, zorder=0, vmin=-10, vmax=20)
+                gdf[wb].plot("winter_stage", ax=ax, zorder=0, norm=norm, cmap=cmap)
             c = waterboards.at[wb, "geometry"].centroid
             ax.text(c.x, c.y, wb.replace(" ", "\n"), ha="center", va="center")
+        nlmod.plot.colorbar_inside(
+            ax=ax,
+            norm=norm,
+            cmap=cmap,
+            bounds=[0.05, 0.55, 0.02, 0.4],
+            label="Summer stage (m NAP)",
+        )
 
 
 @pytest.mark.skip("too slow")