From 142cc222c4a4d3f267107452db4d20bc9fdb3b40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Thu, 8 Jun 2023 08:42:49 +0200 Subject: [PATCH 1/6] Update version.py --- nlmod/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nlmod/version.py b/nlmod/version.py index c7e72605..728ef1b2 100644 --- a/nlmod/version.py +++ b/nlmod/version.py @@ -1,7 +1,7 @@ from importlib import metadata from platform import python_version -__version__ = "0.6.0" +__version__ = "0.6.1b" def show_versions() -> None: From d2432f4ac4494b3ee8581e0d63c5d2d39f542b29 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Wed, 21 Jun 2023 14:35:09 +0200 Subject: [PATCH 2/6] fixes in base.py - reorder coords - add xorigin/yorigin to x,y coords in _get_structured_grid_ds --- nlmod/dims/base.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/nlmod/dims/base.py b/nlmod/dims/base.py index e1569e5e..b4e62e86 100644 --- a/nlmod/dims/base.py +++ b/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( From b75f83abe5128fb8b6f4469b22fb8df8fd43e673 Mon Sep 17 00:00:00 2001 From: dbrakenhoff Date: Wed, 21 Jun 2023 14:35:40 +0200 Subject: [PATCH 3/6] fix setting index in gdf_to_da for vertex grids --- nlmod/dims/grid.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 01b5c4cf..110b5875 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -1062,9 +1062,12 @@ def gdf_to_da( 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] From 18db5c58efffbfd7bbaa51282685a2c447e55b5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ruben=20Calj=C3=A9?= Date: Fri, 23 Jun 2023 20:02:30 +0200 Subject: [PATCH 4/6] Update test_008_waterschappen.py --- tests/test_008_waterschappen.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/tests/test_008_waterschappen.py b/tests/test_008_waterschappen.py index ac1217d0..0ec2b10b 100644 --- a/tests/test_008_waterschappen.py +++ b/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") From e1da3781461bbf2c3d802afb469fd0fbbaf74ec9 Mon Sep 17 00:00:00 2001 From: OnnoEbbens Date: Mon, 26 Jun 2023 13:48:07 +0200 Subject: [PATCH 5/6] fix aggregate nearest for different modelgrids. further deprecating gdf_to_data_array_struc --- docs/examples/06_gridding_vector_data.ipynb | 53 +++++++++---- nlmod/dims/grid.py | 84 ++++++++++++++------- 2 files changed, 94 insertions(+), 43 deletions(-) diff --git a/docs/examples/06_gridding_vector_data.ipynb b/docs/examples/06_gridding_vector_data.ipynb index de58d891..11455287 100644 --- a/docs/examples/06_gridding_vector_data.ipynb +++ b/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,6 +116,7 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ @@ -128,6 +124,7 @@ ] }, { + "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, diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index 110b5875..ccbf7019 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -56,7 +56,8 @@ def xy_to_icell2d(xy, ds): number of the icell2d value of a cell containing the xy point. """ - icell2d = (np.abs(ds.x.data - xy[0]) + np.abs(ds.y.data - xy[1])).argmin().item() + icell2d = (np.abs(ds.x.data - xy[0]) + + np.abs(ds.y.data - xy[1])).argmin().item() return icell2d @@ -354,7 +355,8 @@ def refine( fname, geom_type, level = refinement_feature if not model_coordinates and ds_has_rotation: raise ( - Exception("Converting files to model coordinates not supported") + Exception( + "Converting files to model coordinates not supported") ) g.add_refinement_features(fname, geom_type, level, layers=[0]) elif len(refinement_feature) == 2: @@ -376,8 +378,10 @@ def refine( ) mask = geom_types == geom_type # features = [gdf[mask].unary_union] - features = list(gdf[mask].geometry.explode(index_parts=True)) - g.add_refinement_features(features, geom_type, level, layers=[0]) + features = list( + gdf[mask].geometry.explode(index_parts=True)) + g.add_refinement_features( + features, geom_type, level, layers=[0]) g.build() gridprops = g.get_gridprops_disv() gridprops["area"] = g.get_area() @@ -419,13 +423,15 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", nodata=-1): y = xr.DataArray(xyi[:, 1], dims=("icell2d",)) if method in ["nearest", "linear"]: # resample the entire dataset in one line - ds_out = ds_in.interp(x=x, y=y, method=method, kwargs={"fill_value": None}) + ds_out = ds_in.interp(x=x, y=y, method=method, + kwargs={"fill_value": None}) else: ds_out = xr.Dataset(coords={"layer": ds_in.layer.data, "x": x, "y": y}) # add other variables for data_var in ds_in.data_vars: - data_arr = structured_da_to_ds(ds_in[data_var], ds_out, method=method) + data_arr = structured_da_to_ds( + ds_in[data_var], ds_out, method=method) ds_out[data_var] = data_arr if "area" in gridprops: @@ -523,7 +529,8 @@ def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs): else: x = ds.x y = ds.y - layer_ds = layer_ds.interp(x=x, y=y, method=method, kwargs={"fill_value": None}) + layer_ds = layer_ds.interp( + x=x, y=y, method=method, kwargs={"fill_value": None}) for var in layer_ds.data_vars: ds[var] = layer_ds[var] else: @@ -586,7 +593,8 @@ def col_to_list(col_in, ds, cellids): # 2d vertex grid col_lst = col_in.data[cellids[0]] else: - raise ValueError(f"could not create a column list for col_in={col_in}") + raise ValueError( + f"could not create a column list for col_in={col_in}") else: col_lst = [col_in] * len(cellids[0]) @@ -677,7 +685,8 @@ def lrc_to_reclist( for i_aux in aux: if isinstance(i_aux, str): if "layer" in ds[i_aux].dims and len(cellids) != 3: - cols.append(col_to_list(i_aux, ds, (np.array(layers),) + cellids)) + cols.append(col_to_list( + i_aux, ds, (np.array(layers),) + cellids)) else: cols.append(col_to_list(i_aux, ds, cellids)) else: @@ -772,7 +781,8 @@ def lcid_to_reclist( for i_aux in aux: if isinstance(i_aux, str): if "layer" in ds[i_aux].dims and len(cellids) != 2: - cols.append(col_to_list(i_aux, ds, (np.array(layers),) + cellids)) + cols.append(col_to_list( + i_aux, ds, (np.array(layers),) + cellids)) else: cols.append(col_to_list(i_aux, ds, cellids)) else: @@ -984,7 +994,8 @@ def gdf_to_data_array_struc( # interpolate data if interp_method is not None: - arr = interpolate_gdf_to_array(gdf, gwf, field=field, method=interp_method) + arr = interpolate_gdf_to_array( + gdf, gwf, field=field, method=interp_method) da.values = arr return da @@ -997,7 +1008,8 @@ def gdf_to_data_array_struc( raise ValueError( "multiple geometries in one cell please define aggregation method" ) - gdf_agg = aggregate_vector_per_cell(gdf_cellid, {field: agg_method}, gwf) + gdf_agg = aggregate_vector_per_cell( + gdf_cellid, {field: agg_method}, gwf) else: # aggregation not neccesary gdf_agg = gdf_cellid[[field]] @@ -1058,7 +1070,13 @@ 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]] @@ -1098,7 +1116,8 @@ def interpolate_gdf_to_array(gdf, gwf, field="values", method="nearest"): # check geometry geom_types = gdf.geometry.type.unique() if geom_types[0] != "Point": - raise NotImplementedError("can only use interpolation with point geometries") + raise NotImplementedError( + "can only use interpolation with point geometries") # check field if field not in gdf.columns: @@ -1134,21 +1153,29 @@ def _agg_max_length(gdf, col): def _agg_length_weighted(gdf, col): nanmask = gdf[col].isna() - aw = (gdf.length * gdf[col]).sum(skipna=True) / gdf.loc[~nanmask].length.sum() + aw = (gdf.length * gdf[col]).sum(skipna=True) / \ + gdf.loc[~nanmask].length.sum() 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 @@ -1163,7 +1190,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 @@ -1180,7 +1207,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 @@ -1192,7 +1219,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 @@ -1230,7 +1257,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 @@ -1570,7 +1597,8 @@ def mask_model_edge(ds, idomain): """ # add constant head cells at model boundaries if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: - raise NotImplementedError("model edge not yet calculated for rotated grids") + raise NotImplementedError( + "model edge not yet calculated for rotated grids") # get mask with grid edges xmin = ds["x"] == ds["x"].min() From 405e0ecaaaf2ab1f747cdc9145cf1e7476cafde2 Mon Sep 17 00:00:00 2001 From: OnnoEbbens Date: Mon, 26 Jun 2023 14:30:01 +0200 Subject: [PATCH 6/6] black and version bump --- nlmod/dims/grid.py | 67 +++++++++++++++++----------------------------- nlmod/version.py | 2 +- 2 files changed, 26 insertions(+), 43 deletions(-) diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py index ccbf7019..c74294c4 100644 --- a/nlmod/dims/grid.py +++ b/nlmod/dims/grid.py @@ -56,8 +56,7 @@ def xy_to_icell2d(xy, ds): number of the icell2d value of a cell containing the xy point. """ - icell2d = (np.abs(ds.x.data - xy[0]) + - np.abs(ds.y.data - xy[1])).argmin().item() + icell2d = (np.abs(ds.x.data - xy[0]) + np.abs(ds.y.data - xy[1])).argmin().item() return icell2d @@ -355,8 +354,7 @@ def refine( fname, geom_type, level = refinement_feature if not model_coordinates and ds_has_rotation: raise ( - Exception( - "Converting files to model coordinates not supported") + Exception("Converting files to model coordinates not supported") ) g.add_refinement_features(fname, geom_type, level, layers=[0]) elif len(refinement_feature) == 2: @@ -378,10 +376,8 @@ def refine( ) mask = geom_types == geom_type # features = [gdf[mask].unary_union] - features = list( - gdf[mask].geometry.explode(index_parts=True)) - g.add_refinement_features( - features, geom_type, level, layers=[0]) + features = list(gdf[mask].geometry.explode(index_parts=True)) + g.add_refinement_features(features, geom_type, level, layers=[0]) g.build() gridprops = g.get_gridprops_disv() gridprops["area"] = g.get_area() @@ -423,15 +419,13 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", nodata=-1): y = xr.DataArray(xyi[:, 1], dims=("icell2d",)) if method in ["nearest", "linear"]: # resample the entire dataset in one line - ds_out = ds_in.interp(x=x, y=y, method=method, - kwargs={"fill_value": None}) + ds_out = ds_in.interp(x=x, y=y, method=method, kwargs={"fill_value": None}) else: ds_out = xr.Dataset(coords={"layer": ds_in.layer.data, "x": x, "y": y}) # add other variables for data_var in ds_in.data_vars: - data_arr = structured_da_to_ds( - ds_in[data_var], ds_out, method=method) + data_arr = structured_da_to_ds(ds_in[data_var], ds_out, method=method) ds_out[data_var] = data_arr if "area" in gridprops: @@ -529,8 +523,7 @@ def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs): else: x = ds.x y = ds.y - layer_ds = layer_ds.interp( - x=x, y=y, method=method, kwargs={"fill_value": None}) + layer_ds = layer_ds.interp(x=x, y=y, method=method, kwargs={"fill_value": None}) for var in layer_ds.data_vars: ds[var] = layer_ds[var] else: @@ -593,8 +586,7 @@ def col_to_list(col_in, ds, cellids): # 2d vertex grid col_lst = col_in.data[cellids[0]] else: - raise ValueError( - f"could not create a column list for col_in={col_in}") + raise ValueError(f"could not create a column list for col_in={col_in}") else: col_lst = [col_in] * len(cellids[0]) @@ -685,8 +677,7 @@ def lrc_to_reclist( for i_aux in aux: if isinstance(i_aux, str): if "layer" in ds[i_aux].dims and len(cellids) != 3: - cols.append(col_to_list( - i_aux, ds, (np.array(layers),) + cellids)) + cols.append(col_to_list(i_aux, ds, (np.array(layers),) + cellids)) else: cols.append(col_to_list(i_aux, ds, cellids)) else: @@ -781,8 +772,7 @@ def lcid_to_reclist( for i_aux in aux: if isinstance(i_aux, str): if "layer" in ds[i_aux].dims and len(cellids) != 2: - cols.append(col_to_list( - i_aux, ds, (np.array(layers),) + cellids)) + cols.append(col_to_list(i_aux, ds, (np.array(layers),) + cellids)) else: cols.append(col_to_list(i_aux, ds, cellids)) else: @@ -994,8 +984,7 @@ def gdf_to_data_array_struc( # interpolate data if interp_method is not None: - arr = interpolate_gdf_to_array( - gdf, gwf, field=field, method=interp_method) + arr = interpolate_gdf_to_array(gdf, gwf, field=field, method=interp_method) da.values = arr return da @@ -1008,8 +997,7 @@ def gdf_to_data_array_struc( raise ValueError( "multiple geometries in one cell please define aggregation method" ) - gdf_agg = aggregate_vector_per_cell( - gdf_cellid, {field: agg_method}, gwf) + gdf_agg = aggregate_vector_per_cell(gdf_cellid, {field: agg_method}, gwf) else: # aggregation not neccesary gdf_agg = gdf_cellid[[field]] @@ -1070,13 +1058,13 @@ def gdf_to_da( raise ValueError( "multiple geometries in one cell please define aggregation method" ) - if agg_method in ['nearest']: + 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}) + 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]] @@ -1116,8 +1104,7 @@ def interpolate_gdf_to_array(gdf, gwf, field="values", method="nearest"): # check geometry geom_types = gdf.geometry.type.unique() if geom_types[0] != "Point": - raise NotImplementedError( - "can only use interpolation with point geometries") + raise NotImplementedError("can only use interpolation with point geometries") # check field if field not in gdf.columns: @@ -1153,23 +1140,20 @@ def _agg_max_length(gdf, col): def _agg_length_weighted(gdf, col): nanmask = gdf[col].isna() - aw = (gdf.length * gdf[col]).sum(skipna=True) / \ - gdf.loc[~nanmask].length.sum() + aw = (gdf.length * gdf[col]).sum(skipna=True) / gdf.loc[~nanmask].length.sum() return aw def _agg_nearest(gdf, col, modelgrid): - if modelgrid.grid_type == 'structured': + if modelgrid.grid_type == "structured": cid = gdf["cellid"].values[0] cellcenter = Point( - modelgrid.xcellcenters[0, cid[1]], - modelgrid.ycellcenters[cid[0], 0]) + modelgrid.xcellcenters[0, cid[1]], modelgrid.ycellcenters[cid[0], 0] + ) val = gdf.iloc[gdf.distance(cellcenter).argmin()].loc[col] - elif modelgrid.grid_type == 'vertex': + elif modelgrid.grid_type == "vertex": cid = gdf["cellid"].values[0] - cellcenter = Point( - modelgrid.xcellcenters[cid], - modelgrid.ycellcenters[cid]) + cellcenter = Point(modelgrid.xcellcenters[cid], modelgrid.ycellcenters[cid]) val = gdf.iloc[gdf.distance(cellcenter).argmin()].loc[col] return val @@ -1597,8 +1581,7 @@ def mask_model_edge(ds, idomain): """ # add constant head cells at model boundaries if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0: - raise NotImplementedError( - "model edge not yet calculated for rotated grids") + raise NotImplementedError("model edge not yet calculated for rotated grids") # get mask with grid edges xmin = ds["x"] == ds["x"].min() diff --git a/nlmod/version.py b/nlmod/version.py index 728ef1b2..86efe7c1 100644 --- a/nlmod/version.py +++ b/nlmod/version.py @@ -1,7 +1,7 @@ from importlib import metadata from platform import python_version -__version__ = "0.6.1b" +__version__ = "0.6.1" def show_versions() -> None: