diff --git a/docs/examples/09_schoonhoven.ipynb b/docs/examples/09_schoonhoven.ipynb index 88b723dc..36243bc8 100644 --- a/docs/examples/09_schoonhoven.ipynb +++ b/docs/examples/09_schoonhoven.ipynb @@ -31,9 +31,13 @@ "from nlmod.visualise.netcdf import DatasetCrossSection\n", "from rasterstats import zonal_stats\n", "from shapely.geometry import LineString, Point\n", + "import warnings\n", + "from shapely.errors import ShapelyDeprecationWarning\n", + "\n", + "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n", "\n", "# set the log-level to INFO, so more information is shown (compared to the default setting of WARNING)\n", - "logging.basicConfig(level=logging.INFO)\n" + "logging.basicConfig(level=logging.INFO)" ] }, { @@ -157,7 +161,7 @@ "metadata": {}, "outputs": [], "source": [ - "pg = nlmod.gwf.surface_water.download_level_areas(bgt, extent=extent)\n" + "pg = nlmod.gwf.surface_water.download_level_areas(bgt, extent=extent)" ] }, { @@ -179,7 +183,7 @@ "f, ax = nlmod.plot.get_map(extent)\n", "bgt.plot(color=\"k\", ax=ax)\n", "for wb in pg:\n", - " pg[wb].plot(\"summer_stage\", ax=ax, vmin=-3, vmax=1, zorder=0)\n" + " pg[wb].plot(\"summer_stage\", ax=ax, vmin=-3, vmax=1, zorder=0)" ] }, { @@ -221,7 +225,7 @@ "norm = matplotlib.colors.Normalize(vmin=-3, vmax=1)\n", "cmap = \"viridis\"\n", "bgt.plot(\"summer_stage\", ax=ax, norm=norm, cmap=cmap)\n", - "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)\n" + "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)" ] }, { @@ -241,7 +245,7 @@ "source": [ "f, ax = nlmod.plot.get_map(extent)\n", "bgt.plot(\"ahn_min\", ax=ax, norm=norm, cmap=cmap)\n", - "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)\n" + "nlmod.plot.colorbar_inside(norm=norm, cmap=cmap)" ] }, { @@ -261,7 +265,7 @@ "outputs": [], "source": [ "regis = nlmod.read.get_regis(extent, cachedir=cachedir, cachename=\"regis.nc\")\n", - "regis\n" + "regis" ] }, { @@ -280,7 +284,7 @@ "outputs": [], "source": [ "ds = nlmod.to_model_ds(regis, model_name, model_ws, delr=100.0, delc=100.0)\n", - "ds\n" + "ds" ] }, { @@ -416,7 +420,7 @@ "bgt_grid[\"cond\"] = bgt_grid.area / bed_resistance\n", "mask = bgt_grid[\"bronhouder\"] == \"L0002\"\n", "lek = bgt_grid[mask]\n", - "bgt_grid = bgt_grid[~mask]\n" + "bgt_grid = bgt_grid[~mask]" ] }, { @@ -520,7 +524,7 @@ "cbar = nlmod.plot.colorbar_inside(pc)\n", "for label in cbar.ax.yaxis.get_ticklabels():\n", " label.set_bbox(dict(facecolor=\"w\", alpha=0.5))\n", - "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")\n" + "bgt.plot(ax=ax, edgecolor=\"k\", facecolor=\"none\")" ] }, { @@ -550,7 +554,7 @@ " label.set_bbox(dict(facecolor=\"w\", alpha=0.5))\n", "dcs.plot_grid()\n", "dcs.plot_layers(alpha=0.0, min_label_area=1000)\n", - "f.tight_layout(pad=0.0)\n" + "f.tight_layout(pad=0.0)" ] }, { @@ -577,7 +581,7 @@ " head_point = head.interp(x=x, y=y, method=\"nearest\")\n", "# only keep layers that are active at this location\n", "head_point = head_point[:, ~head_point.isnull().all(\"time\")]\n", - "head_point.plot.line(hue=\"layer\", size=(10));\n" + "head_point.plot.line(hue=\"layer\", size=10);" ] }, { @@ -617,7 +621,7 @@ " ha=\"center\",\n", " va=\"top\",\n", " transform=ax.transAxes,\n", - " )\n" + " )" ] }, { @@ -687,7 +691,7 @@ " pf = pdata[pdata['particleid']==pid]\n", " ax.plot(pf['x'],pf['y'], color=\"k\", linewidth=0.5)\n", "ax.plot(pf['x'],pf['y'], color=\"k\", linewidth=0.5, label='pathline')\n", - "bgt.plot(ax=ax, edgecolor=\"blue\", facecolor=\"none\");\n" + "bgt.plot(ax=ax, edgecolor=\"blue\", facecolor=\"none\");" ] }, { @@ -713,16 +717,8 @@ "dcs.plot_grid()\n", "# add labels with layer names\n", "dcs.plot_layers(alpha=0.0, min_label_area=1000)\n", - "f.tight_layout(pad=0.0)\n" + "f.tight_layout(pad=0.0)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "008d54c3-3ef7-4766-8447-b7837943db6e", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/docs/examples/10_modpath.ipynb b/docs/examples/10_modpath.ipynb index 343b4a1b..6591102d 100644 --- a/docs/examples/10_modpath.ipynb +++ b/docs/examples/10_modpath.ipynb @@ -92,7 +92,7 @@ "sim = flopy.mf6.MFSimulation.load(\n", " \"mfsim.nam\", sim_ws=model_ws, exe_name=exe_name\n", ")\n", - "gwf = sim.get_model(model_name=model_name)\n" + "gwf = sim.get_model(model_name=model_name)" ] }, { @@ -131,7 +131,7 @@ "pg = nlmod.modpath.pg_from_fdt(nodes)\n", "\n", "# create the modpath simulation file\n", - "mpsim = nlmod.modpath.sim(mpf, pg, \"backward\", gwf=gwf)\n" + "mpsim = nlmod.modpath.sim(mpf, pg, \"backward\", gwf=gwf)" ] }, { @@ -178,7 +178,7 @@ " color=\"red\",\n", ")\n", "ax.set_title(f\"pathlines\")\n", - "ax.legend()\n" + "ax.legend()" ] }, { @@ -200,7 +200,7 @@ "ax.set_ylabel(\"distance [m]\")\n", "ax.set_xlabel(\"time [year]\")\n", "ax.set_title(\"distance travelled per particle\")\n", - "ax.grid()\n" + "ax.grid()" ] }, { @@ -233,7 +233,7 @@ "pg = nlmod.modpath.pg_from_fdt(nodes)\n", "\n", "# create the modpath simulation file\n", - "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\")\n" + "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\")" ] }, { @@ -288,7 +288,7 @@ " ax.set_ylim(498700, 499300)\n", " elif i == 2:\n", " ax.set_xlim(101200, 101700)\n", - " ax.set_ylim(496300, 496700)\n" + " ax.set_ylim(496300, 496700)" ] }, { @@ -310,7 +310,7 @@ "ax.set_ylabel(\"distance [m]\")\n", "ax.set_xlabel(\"time [year]\")\n", "ax.set_title(\"distance travelled per particle\")\n", - "ax.grid()\n" + "ax.grid()" ] }, { @@ -348,7 +348,7 @@ "pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)\n", "\n", "# create the modpath simulation file\n", - "mpsim = nlmod.modpath.sim(mpf, pg, \"backward\", gwf=gwf)\n" + "mpsim = nlmod.modpath.sim(mpf, pg, \"backward\", gwf=gwf)" ] }, { @@ -404,7 +404,7 @@ " ax.set_ylim(498300, 499300)\n", " elif i == 2:\n", " ax.set_xlim(101000, 102000)\n", - " ax.set_ylim(496300, 497300)\n" + " ax.set_ylim(496300, 497300)" ] }, { @@ -427,7 +427,7 @@ "ax.set_ylabel(\"distance [m]\")\n", "ax.set_xlabel(\"time [year]\")\n", "ax.set_title(\"distance travelled per particle\")\n", - "ax.grid()\n" + "ax.grid()" ] }, { @@ -458,7 +458,7 @@ "pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)\n", "\n", "# create the modpath simulation file\n", - "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf, stoptime=10 * 365)\n" + "mpsim = nlmod.modpath.sim(mpf, pg, \"forward\", gwf=gwf, stoptime=10 * 365)" ] }, { @@ -513,7 +513,7 @@ " ax.set_ylim(498300, 499300)\n", " elif i == 2:\n", " ax.set_xlim(101000, 102000)\n", - " ax.set_ylim(496300, 497300)\n" + " ax.set_ylim(496300, 497300)" ] }, { @@ -536,7 +536,7 @@ "ax.set_ylabel(\"distance [m]\")\n", "ax.set_xlabel(\"time [year]\")\n", "ax.set_title(\"distance travelled per particle\")\n", - "ax.grid()\n" + "ax.grid()" ] } ], diff --git a/docs/examples/13_plot_methods.ipynb b/docs/examples/13_plot_methods.ipynb new file mode 100644 index 00000000..87db3c6b --- /dev/null +++ b/docs/examples/13_plot_methods.ipynb @@ -0,0 +1,228 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a3f62867", + "metadata": {}, + "source": [ + "# Plot methods\n", + "This notebook shows the plot methods that are available in nlmod. Most plot methods use a model Dataset as input, which is an xarray Dataset with some required variables and attributes. There are some plot methods in flopy as well, whcih require a grounwwater flow model, or a modelgrid?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9ae5a49", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import logging\n", + "import matplotlib.pyplot as plt\n", + "import xarray as xr\n", + "import flopy\n", + "import nlmod\n", + "from nlmod.visualise.netcdf import DatasetCrossSection\n", + "\n", + "# set the log-level to INFO, so more information is shown (compared to the default setting of WARNING)\n", + "logging.basicConfig(level=logging.INFO)" + ] + }, + { + "cell_type": "markdown", + "id": "f5ca2bf5", + "metadata": {}, + "source": [ + "First we read a fully run model, from notebook 9. Please run this notebook first." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ca958d0", + "metadata": {}, + "outputs": [], + "source": [ + "model_name = \"Schoonhoven\"\n", + "model_ws = \"model9\"\n", + "figdir, cachedir = nlmod.util.get_model_dirs(model_ws)\n", + "ds = xr.open_dataset(os.path.join(cachedir, 'full_ds.nc'))\n", + "# add calculated heads\n", + "ds['head'] = nlmod.util.get_heads_dataarray(ds)\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "id": "3909cdd5", + "metadata": {}, + "source": [ + "For the flopy plot-methods we need a modelgrid object. We generate this from the model Dataset using the method. nlmod.mgrid.modelgrid_from_ds()." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a0b7349", + "metadata": {}, + "outputs": [], + "source": [ + "modelgrid = nlmod.mgrid.modelgrid_from_ds(ds)\n", + "modelgrid" + ] + }, + { + "cell_type": "markdown", + "id": "f594a18d", + "metadata": {}, + "source": [ + "## Maps\n", + "We can plot variables on a map using nlmod.plot.data_array(). We can also use the PlotMapView-class from flopy, and plot an array using the plot_array method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "909de9e4", + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = nlmod.plot.get_map(ds.extent, ncols=2)\n", + "\n", + "# plot using nlmod\n", + "pc = nlmod.plot.data_array(ds['top'], ds=ds, ax=ax[0])\n", + "\n", + "# plot using flopy\n", + "pmv = flopy.plot.PlotMapView(modelgrid=modelgrid, ax=ax[1])\n", + "pmv.plot_array(ds['top'])" + ] + }, + { + "cell_type": "markdown", + "id": "98368c1d", + "metadata": {}, + "source": [ + "## Cross-sections\n", + "We can also plot cross-sections, either with DatasetCrossSection in nlmod, or using the PlotCrossSection class of flopy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32a89356", + "metadata": {}, + "outputs": [], + "source": [ + "y = (ds.extent[2] + ds.extent[3]) / 2 + 0.1\n", + "line = [(ds.extent[0], y), (ds.extent[1], y)]\n", + "zmin = -100.\n", + "zmax = 10." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "044aa8f5", + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = plt.subplots(figsize=(10,5), nrows=2)\n", + "\n", + "# plot using nlmod\n", + "dcs = DatasetCrossSection(ds, line=line, zmin=zmin, zmax=zmax, ax=ax[0])\n", + "dcs.plot_array(ds['kh'])\n", + "\n", + "# plot using flopy\n", + "pcs = flopy.plot.PlotCrossSection(modelgrid=modelgrid, line={'line':line}, ax=ax[1])\n", + "pcs.plot_array(ds['kh'])\n", + "pcs.ax.set_ylim((zmin, zmax))" + ] + }, + { + "cell_type": "markdown", + "id": "b1f49969", + "metadata": {}, + "source": [ + "## Time series\n", + "For time series we use the functionality of xarray, as we have read the heads in a xarray DataArray." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ee45ab1", + "metadata": {}, + "outputs": [], + "source": [ + "x = 118228\n", + "y = 439870\n", + "cellid = modelgrid.intersect(x=x, y=y)\n", + "if isinstance(cellid, int):\n", + " head_point = ds['head'].loc[:, :, cellid]\n", + "else:\n", + " head_point = ds['head'].loc[:, :, cellid[0], cellid[1]]\n", + "# only keep layers that are active at this location\n", + "head_point = head_point[:, ~head_point.isnull().all(\"time\")]\n", + "head_point.plot.line(hue=\"layer\", size=10);" + ] + }, + { + "cell_type": "markdown", + "id": "d56cc81c", + "metadata": {}, + "source": [ + "We can also use pandas to plot the heads. First transform the data to a Pandas DataFrame." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b129675", + "metadata": {}, + "outputs": [], + "source": [ + "df = head_point.to_pandas()\n", + "df" + ] + }, + { + "cell_type": "markdown", + "id": "02112ab3", + "metadata": {}, + "source": [ + "And then plot this DataFrame." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd9d1f42", + "metadata": {}, + "outputs": [], + "source": [ + "df.plot(figsize=(10,10))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nlmod/mdims/mgrid.py b/nlmod/mdims/mgrid.py index f3df835f..483bb80d 100644 --- a/nlmod/mdims/mgrid.py +++ b/nlmod/mdims/mgrid.py @@ -63,7 +63,7 @@ def xy_to_icell2d(xy, ds): return icell2d -def modelgrid_from_ds(ds, rotated=True, **kwargs): +def modelgrid_from_ds(ds, rotated=True, nlay=None, top=None, botm=None, **kwargs): """Get flopy modelgrid from ds. Parameters @@ -88,7 +88,22 @@ def modelgrid_from_ds(ds, rotated=True, **kwargs): xoff = 0.0 yoff = 0.0 angrot = 0.0 - + if top is None and "top" in ds: + top = ds["top"].data + if botm is None and "botm" in ds: + botm = ds["botm"].data + if nlay is None: + if "layer" in ds: + nlay = len(ds.layer) + elif botm is not None: + nlay = len(botm) + + if nlay is not None and botm is not None and nlay < len(botm): + botm = botm[:nlay] + + kwargs = dict( + xoff=xoff, yoff=yoff, angrot=angrot, nlay=nlay, top=top, botm=botm, **kwargs + ) if ds.gridtype == "structured": if not isinstance(ds.extent, (tuple, list, np.ndarray)): raise TypeError( @@ -99,9 +114,6 @@ def modelgrid_from_ds(ds, rotated=True, **kwargs): modelgrid = StructuredGrid( delc=delc, delr=delr, - xoff=xoff, - yoff=yoff, - angrot=angrot, **kwargs, ) elif ds.gridtype == "vertex": @@ -110,9 +122,6 @@ def modelgrid_from_ds(ds, rotated=True, **kwargs): modelgrid = VertexGrid( vertices=vertices, cell2d=cell2d, - xoff=xoff, - yoff=yoff, - angrot=angrot, **kwargs, ) return modelgrid @@ -168,11 +177,17 @@ def get_cell2d_from_ds(ds): x = ds["x"].data y = ds["y"].data icvert = ds["icvert"].data + if "_FillValue" in ds["icvert"].attrs: + nodata = ds["icvert"].attrs["_FillValue"] + else: + nodata = -1 + icvert = icvert.copy() + icvert[np.isnan(icvert)] = nodata + icvert = icvert.astype(int) cell2d = [] - nodata = ds["icvert"].attrs["_FillValue"] for i, cid in enumerate(icell2d): - mask = ds["icvert"].data[i] != nodata - cell2d.append((cid, x[i], y[i], mask.sum(), *icvert[i][mask])) + mask = icvert[i] != nodata + cell2d.append((cid, x[i], y[i], mask.sum(), *icvert[i, mask])) return cell2d @@ -238,21 +253,13 @@ def refine( g = Gridgen(dis, model_ws=model_ws, exe_name=exe_name) else: # create a modelgrid with only one layer, to speed up Gridgen - top = ds["top"].values - botm = ds["botm"].values[[0]] - modelgrid = modelgrid_from_ds( - ds, rotated=False, nlay=1, top=top, botm=botm - ) + modelgrid = modelgrid_from_ds(ds, rotated=False, nlay=1) g = Gridgen(modelgrid, model_ws=model_ws, exe_name=exe_name) ds_has_rotation = "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0 if model_coordinates: if not ds_has_rotation: - raise ( - Exception( - "The supplied shapes need to be in realworld coordinates" - ) - ) + raise (Exception("The supplied shapes need to be in realworld coordinates")) elif ds_has_rotation: affine_matrix = get_affine_world_to_mod(ds).to_shapely() @@ -263,9 +270,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: @@ -288,9 +293,7 @@ def refine( mask = geom_types == geom_type # features = [gdf[mask].unary_union] features = list(gdf[mask].geometry.explode()) - g.add_refinement_features( - features, geom_type, level, layers=[0] - ) + g.add_refinement_features(features, geom_type, level, layers=[0]) g.build() gridprops = g.get_gridprops_disv() gridprops["area"] = g.get_area() @@ -396,25 +399,20 @@ def col_to_list(col_in, ds, cellids): elif len(cellids) == 2: # 2d grid or vertex 3d grid col_lst = [ - col_in.data[row, col] - for row, col in zip(cellids[0], cellids[1]) + col_in.data[row, col] for row, col in zip(cellids[0], cellids[1]) ] elif len(cellids) == 1: # 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]) return col_lst -def lrc_to_reclist( - layers, rows, columns, cellids, ds, col1=None, col2=None, col3=None -): +def lrc_to_reclist(layers, rows, columns, cellids, ds, col1=None, col2=None, col3=None): """Create a reclist for stress period data from a set of cellids. Used for structured grids. @@ -484,13 +482,9 @@ def lrc_to_reclist( col1_lst = col_to_list(col1, ds, cellids) col2_lst = col_to_list(col2, ds, cellids) col3_lst = col_to_list(col3, ds, cellids) - reclist = list( - zip(zip(layers, rows, columns), col1_lst, col2_lst, col3_lst) - ) + reclist = list(zip(zip(layers, rows, columns), col1_lst, col2_lst, col3_lst)) else: - raise ValueError( - "invalid combination of values for col1, col2 and col3" - ) + raise ValueError("invalid combination of values for col1, col2 and col3") return reclist @@ -560,13 +554,9 @@ def lcid_to_reclist(layers, cellids, ds, col1=None, col2=None, col3=None): col1_lst = col_to_list(col1, ds, cellids) col2_lst = col_to_list(col2, ds, cellids) col3_lst = col_to_list(col3, ds, cellids) - reclist = list( - zip(zip(layers, cellids[-1]), col1_lst, col2_lst, col3_lst) - ) + reclist = list(zip(zip(layers, cellids[-1]), col1_lst, col2_lst, col3_lst)) else: - raise ValueError( - "invalid combination of values for col1, col2 and col3" - ) + raise ValueError("invalid combination of values for col1, col2 and col3") return reclist @@ -653,9 +643,7 @@ def da_to_reclist( layers = cellids[0] rows = cellids[1] columns = cellids[2] - return lrc_to_reclist( - layers, rows, columns, cellids, ds, col1, col2, col3 - ) + return lrc_to_reclist(layers, rows, columns, cellids, ds, col1, col2, col3) else: if first_active_layer: fal = get_first_active_layer(ds) @@ -679,9 +667,7 @@ def da_to_reclist( rows = cellids[-2] columns = cellids[-1] - return lrc_to_reclist( - layers, rows, columns, cellids, ds, col1, col2, col3 - ) + return lrc_to_reclist(layers, rows, columns, cellids, ds, col1, col2, col3) def polygon_to_area(modelgrid, polygon, da, gridtype="structured"): @@ -766,9 +752,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 @@ -781,9 +765,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]] @@ -935,9 +917,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: @@ -973,9 +953,7 @@ 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 @@ -1058,9 +1036,7 @@ def aggregate_vector_per_cell(gdf, fields_methods, gwf=None): if ("Polygon" in geom_types) or ("MultiPolygon" in geom_types): pass else: - raise TypeError( - "can only use area methods with polygon geometries" - ) + raise TypeError("can only use area methods with polygon geometries") # check fields missing_cols = set(fields_methods.keys()).difference(gdf.columns) @@ -1107,9 +1083,7 @@ def gdf_to_bool_data_array(gdf, mfgrid, ds): elif ds.gridtype == "vertex": da = util.get_da_from_da_ds(ds, dims=("icell2d",), data=0) else: - raise ValueError( - "function only support structured or vertex gridtypes" - ) + raise ValueError("function only support structured or vertex gridtypes") if isinstance(gdf, gpd.GeoDataFrame): geoms = gdf.geometry.values @@ -1122,8 +1096,8 @@ def gdf_to_bool_data_array(gdf, mfgrid, ds): ncol = mfgrid.ncol for cid in cids: if version.parse(flopy.__version__) < version.parse("3.3.6"): - i, j = cid - else: + i, j = cid + else: # TODO: temporary fix until flopy intersect on structured # grid returns row, col again. i = int((cid) / ncol) @@ -1256,9 +1230,7 @@ def get_thickness_from_topbot(top, bot): elif bot.ndim == 2: thickness = util.get_da_from_da_ds(bot, dims=("layer", "icell2d")) else: - raise ValueError( - "function only support structured or vertex gridtypes" - ) + raise ValueError("function only support structured or vertex gridtypes") for lay in range(len(bot)): if lay == 0: @@ -1269,9 +1241,7 @@ def get_thickness_from_topbot(top, bot): return thickness -def get_vertices_arr( - ds, modelgrid=None, vert_per_cid=4, epsilon=0, rotated=False -): +def get_vertices_arr(ds, modelgrid=None, vert_per_cid=4, epsilon=0, rotated=False): """get vertices of a vertex modelgrid from a ds or the modelgrid. Only return the 4 corners of each cell and not the corners of adjacent cells thus limiting the vertices per cell to 4 points. @@ -1411,9 +1381,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/visualise/plots.py b/nlmod/visualise/plots.py index 9f3905c4..e78b7c65 100644 --- a/nlmod/visualise/plots.py +++ b/nlmod/visualise/plots.py @@ -416,11 +416,16 @@ def get_patches(ds, rotated=False): affine = get_affine_mod_to_world(ds) xy[:, 0], xy[:, 1] = affine * (xy[:, 0], xy[:, 1]) icvert = ds["icvert"].data - nodata = ds["icvert"].attrs["_FillValue"] - patches = [ - Polygon(xy[icvert[icell2d, icvert[icell2d] != nodata]]) - for icell2d in ds.icell2d.data - ] + if "_FillValue" in ds["icvert"].attrs: + nodata = ds["icvert"].attrs["_FillValue"] + else: + nodata = -1 + icvert = icvert.copy() + icvert[np.isnan(icvert)] = nodata + icvert = icvert.astype(int) + patches = [] + for icell2d in ds.icell2d.data: + patches.append(Polygon(xy[icvert[icell2d, icvert[icell2d] != nodata]])) return patches