diff --git a/autotest/test_binaryfile.py b/autotest/test_binaryfile.py index 209e7c456f..b224f2209b 100644 --- a/autotest/test_binaryfile.py +++ b/autotest/test_binaryfile.py @@ -274,7 +274,7 @@ def test_load_binary_head_file(example_data_path): def test_plot_binary_head_file(example_data_path): hf = HeadFile(example_data_path / "freyberg" / "freyberg.githds") - hf.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) + hf.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) assert isinstance(hf.plot(), Axes) plt.close() diff --git a/autotest/test_formattedfile.py b/autotest/test_formattedfile.py index 0d31ed308c..26ac20f451 100644 --- a/autotest/test_formattedfile.py +++ b/autotest/test_formattedfile.py @@ -69,7 +69,7 @@ def test_headfile_build_index(example_data_path): def test_formattedfile_reference(example_data_path): h = FormattedHeadFile(example_data_path / "mf2005_test" / "test1tr.githds") assert isinstance(h, FormattedHeadFile) - h.mg.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) + h.modelgrid.set_coord_info(xoff=1000.0, yoff=200.0, angrot=15.0) assert isinstance(h.plot(masked_values=[6999.000]), Axes) plt.close() diff --git a/autotest/test_geospatial_util.py b/autotest/test_geospatial_util.py index decbcf3706..135d8b6957 100644 --- a/autotest/test_geospatial_util.py +++ b/autotest/test_geospatial_util.py @@ -374,7 +374,7 @@ def test_polygon_collection(polygon, poly_w_hole, multipolygon): points = gc1.points geojson = gc1.geojson fp_geo = gc1.flopy_geometry - gdf = gc1.geo_dataframe + gdf = gc1.geodataframe collections = [shp, shply, points, geojson, fp_geo, gdf] for col in collections: @@ -408,7 +408,7 @@ def test_point_collection(point, multipoint): points = gc1.points geojson = gc1.geojson fp_geo = gc1.flopy_geometry - gdf = gc1.geo_dataframe + gdf = gc1.geodataframe collections = [shp, shply, points, geojson, fp_geo, gdf] for col in collections: @@ -436,7 +436,7 @@ def test_linestring_collection(linestring, multilinestring): points = gc1.points geojson = gc1.geojson fp_geo = gc1.flopy_geometry - gdf = gc1.geo_dataframe + gdf = gc1.geodataframe collections = [shp, shply, points, geojson, fp_geo, gdf] for col in collections: @@ -481,7 +481,7 @@ def test_mixed_collection( points = gc1.points geojson = gc1.geojson fp_geo = gc1.flopy_geometry - gdf = gc1.geo_dataframe + gdf = gc1.geodataframe collections = [shp, shply, lshply, points, geojson, fp_geo, gdf] for col in collections: @@ -525,7 +525,7 @@ def test_geopandas_dtypes( col = Collection(col) gc1 = GeoSpatialCollection(col) - gdf = gc1.geo_dataframe + gdf = gc1.geodataframe collections = [gdf, gdf.geometry, gdf.geometry.values] for col in collections: diff --git a/autotest/test_grid.py b/autotest/test_grid.py index 6f1211e60b..d5cea3c036 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -1585,14 +1585,14 @@ def test_unstructured_convert(unstructured_grid): @requires_pkg("geopandas") -def test_geo_dataframe(structured_grid, vertex_grid, unstructured_grid): +def test_geodataframe(structured_grid, vertex_grid, unstructured_grid): geopandas = import_optional_dependency("geopandas") grids = (structured_grid, vertex_grid, unstructured_grid) for grid in grids: - gdf = grid.geo_dataframe + gdf = grid.to_geodataframe() if not isinstance(gdf, geopandas.GeoDataFrame): - raise TypeError("geo_dataframe not returning GeoDataFrame object") + raise TypeError("geodataframe not returning GeoDataFrame object") geoms = gdf.geometry.values for node, geom in enumerate(geoms): diff --git a/autotest/test_modpathfile.py b/autotest/test_modpathfile.py index 652ad4f8b5..883d11f1b7 100644 --- a/autotest/test_modpathfile.py +++ b/autotest/test_modpathfile.py @@ -308,9 +308,9 @@ def test_get_destination_endpoint_data( @pytest.mark.parametrize("longfieldname", [True, False]) @requires_exe("mf6", "mp7") -@requires_pkg("pyshp", "shapely", name_map={"pyshp": "shapefile"}) +@requires_pkg("geopandas", "shapely") def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): - from shapefile import Reader + import geopandas as gpd # setup and run model, then copy outputs to function_tmpdir sim, forward_model_name, _, _, _ = mp7_small @@ -330,13 +330,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): # define shapefile path shp_file = ws / "pathlines.shp" - # add a column to the pathline recarray - fieldname = "newfield" + ("longname" if longfieldname else "") - fieldval = "x" - pathlines = [ - rfn.append_fields(pl, fieldname, list(repeat(fieldval, len(pl))), dtypes="|U1") - for pl in pathlines - ] + pline_names = [name[:10] for name in pathlines[0].dtype.names] # write the pathline recarray to shapefile pathline_file.write_shapefile( @@ -350,8 +344,7 @@ def test_write_shapefile(function_tmpdir, mp7_small, longfieldname): assert shp_file.is_file() # load shapefile - with Reader(shp_file) as reader: - fieldnames = [f[0] for f in reader.fields[1:]] - fieldname = "newfiname_" if longfieldname else fieldname - assert fieldname in fieldnames - assert all(r[fieldname] == fieldval for r in reader.iterRecords()) + gdf = gpd.read_file(shp_file) + fieldnames = list(gdf) + for fname in pline_names: + assert fname in fieldnames diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py index 1d59801b05..d450ea6dd0 100644 --- a/flopy/discretization/grid.py +++ b/flopy/discretization/grid.py @@ -593,7 +593,7 @@ def xyzvertices(self): def cross_section_vertices(self): return self.xyzvertices[0], self.xyzvertices[1] - def geo_dataframe(self, features, featuretype="Polygon"): + def to_geodataframe(self, features, featuretype="Polygon"): """ Method returns a geopandas GeoDataFrame of the Grid @@ -606,15 +606,14 @@ def geo_dataframe(self, features, featuretype="Polygon"): gc = GeoSpatialCollection( features, shapetype=[featuretype for _ in range(len(features))] ) - gdf = gc.geo_dataframe + gdf = gc.geodataframe gdf["node"] = gdf.index + 1 if self.crs is not None: gdf = gdf.set_crs(crs=self.crs) return gdf - @property - def grid_line_geo_dataframe(self): + def grid_line_geodataframe(self): """ Method to get a GeoDataFrame of grid lines @@ -622,7 +621,7 @@ def grid_line_geo_dataframe(self): ------- GeoDataFrame """ - gdf = self.geo_dataframe(self.grid_lines, featuretype="LineString") + gdf = self.to_geodataframe(self.grid_lines, featuretype="LineString") gdf = gdf.rename(columns={"node": "number"}) return gdf diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 0d2df58fdf..3b00ee10ab 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -759,8 +759,7 @@ def map_polygons(self): return self._polygons - @property - def geo_dataframe(self): + def to_geodataframe(self): """ Returns a geopandas GeoDataFrame of the model grid @@ -769,11 +768,28 @@ def geo_dataframe(self): GeoDataFrame """ polys = [[list(zip(*i))] for i in zip(*self.cross_section_vertices)] - gdf = super().geo_dataframe(polys) + gdf = super().to_geodataframe(polys) gdf["row"] = sorted(list(range(1, self.nrow + 1)) * self.ncol) gdf["col"] = list(range(1, self.ncol + 1)) * self.nrow return gdf + @property + def geo_dataframe(self): + """ + Returns a geopandas GeoDataFrame of the model grid + + Returns + ------- + GeoDataFrame + """ + import warnings + + warnings.warn( + "geo_dataframe has been deprecated, use to_geodataframe() instead", + DeprecationWarning, + ) + return self.to_geodataframe() + def convert_grid(self, factor): """ Method to scale the model grid based on user supplied scale factors diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py index f34764dd6f..2f46f7dc1e 100644 --- a/flopy/discretization/unstructuredgrid.py +++ b/flopy/discretization/unstructuredgrid.py @@ -580,8 +580,7 @@ def map_polygons(self): return copy.copy(self._polygons) - @property - def geo_dataframe(self): + def to_geodataframe(self): """ Returns a geopandas GeoDataFrame of the model grid @@ -590,9 +589,26 @@ def geo_dataframe(self): GeoDataFrame """ polys = [[self.get_cell_vertices(nn)] for nn in range(self.nnodes)] - gdf = super().geo_dataframe(polys) + gdf = super().to_geodataframe(polys) return gdf + @property + def geo_dataframe(self): + """ + Returns a geopandas GeoDataFrame of the model grid + + Returns + ------- + GeoDataFrame + """ + import warnings + + warnings.warn( + "geo_dataframe has been deprecated, use to_geodataframe() instead", + DeprecationWarning, + ) + return self.to_geodataframe() + def neighbors(self, node=None, **kwargs): """ Method to get nearest neighbors of a cell diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py index dc5be87109..1cb79a2c0d 100644 --- a/flopy/discretization/vertexgrid.py +++ b/flopy/discretization/vertexgrid.py @@ -300,8 +300,7 @@ def map_polygons(self): return copy.copy(self._polygons) - @property - def geo_dataframe(self): + def to_geodataframe(self): """ Returns a geopandas GeoDataFrame of the model grid @@ -313,9 +312,26 @@ def geo_dataframe(self): featuretype = "Polygon" if self._cell1d is not None: featuretype = "multilinestring" - gdf = super().geo_dataframe(cells, featuretype) + gdf = super().to_geodataframe(cells, featuretype) return gdf + @property + def geo_dataframe(self): + """ + Returns a geopandas GeoDataFrame of the model grid + + Returns + ------- + GeoDataFrame + """ + import warnings + + warnings.warn( + "geo_dataframe has been deprecated, use to_geodataframe() instead", + DeprecationWarning, + ) + return self.to_geodataframe() + def convert_grid(self, factor): """ Method to scale the model grid based on user supplied scale factors diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py index 3c0e188279..a553ba5a8c 100644 --- a/flopy/export/shapefile_utils.py +++ b/flopy/export/shapefile_utils.py @@ -17,12 +17,11 @@ import numpy as np from ..datbase import DataInterface, DataType -from ..discretization.grid import Grid from ..utils import Util3d, flopy_io, import_optional_dependency from ..utils.crs import get_crs -def write_gridlines_shapefile(filename: Union[str, PathLike], mg): +def write_gridlines_shapefile(filename: Union[str, PathLike], modelgrid): """ Write a polyline shapefile of the grid lines - a lightweight alternative to polygons. @@ -31,7 +30,7 @@ def write_gridlines_shapefile(filename: Union[str, PathLike], mg): ---------- filename : str or PathLike path of the shapefile to write - mg : flopy.discretization.grid.Grid object + modelgrid : flopy.discretization.grid.Grid object flopy model grid Returns @@ -39,29 +38,26 @@ def write_gridlines_shapefile(filename: Union[str, PathLike], mg): None """ - shapefile = import_optional_dependency("shapefile") - if not isinstance(mg, Grid): + warnings.warn( + "the write_gridlines_shapefile() utility will be deprecated and is " + "replaced by modelgrid.to_geodataframe()", + DeprecationWarning, + ) + + from ..discretization.grid import Grid + + if not isinstance(modelgrid, Grid): raise ValueError( - f"'mg' must be a flopy Grid subclass instance; found '{type(mg)}'" + f"'modelgrid' must be a flopy Grid subclass instance; found '{type(modelgrid)}'" ) - wr = shapefile.Writer(str(filename), shapeType=shapefile.POLYLINE) - wr.field("number", "N", 18, 0) - grid_lines = mg.grid_lines - for i, line in enumerate(grid_lines): - wr.line([line]) - wr.record(i) - wr.close() - try: - write_prj(filename, modelgrid=mg) - except ImportError: - pass - return + gdf = modelgrid.grid_line_geodataframe() + gdf.write_file(filename) def write_grid_shapefile( - path: Union[str, PathLike], - mg, + filename: Union[str, PathLike], + modelgrid, array_dict, nan_val=np.nan, crs=None, @@ -74,9 +70,9 @@ def write_grid_shapefile( Parameters ---------- - path : str or PathLike + filename : str or PathLike shapefile file path - mg : flopy.discretization.grid.Grid object + modelgrid : flopy.discretization.grid.Grid object flopy model grid array_dict : dict dictionary of model input arrays @@ -105,127 +101,53 @@ def write_grid_shapefile( None """ - shapefile = import_optional_dependency("shapefile") - w = shapefile.Writer(str(path), shapeType=shapefile.POLYGON) - w.autoBalance = 1 + warnings.warn( + "the write_grid_shapefile() utility will be deprecated and is " + "replaced by modelgrid.to_geodataframe()", + DeprecationWarning, + ) + + from ..discretization.grid import Grid - if not isinstance(mg, Grid): + if not isinstance(modelgrid, Grid): raise ValueError( - f"'mg' must be a flopy Grid subclass instance; found '{type(mg)}'" + f"'modelgrid' must be a flopy Grid subclass instance; found '{type(modelgrid)}'" ) - elif mg.grid_type == "structured": - verts = [ - mg.get_cell_vertices(i, j) for i in range(mg.nrow) for j in range(mg.ncol) - ] - elif mg.grid_type == "vertex": - verts = [mg.get_cell_vertices(cellid) for cellid in range(mg.ncpl)] - elif mg.grid_type == "unstructured": - verts = [mg.get_cell_vertices(cellid) for cellid in range(mg.nnodes)] - else: - raise NotImplementedError(f"Grid type {mg.grid_type} not supported.") - - # set up the attribute fields and arrays of attributes - if mg.grid_type == "structured": - names = ["node", "row", "column"] + list(array_dict.keys()) - dtypes = [ - ("node", np.dtype("int")), - ("row", np.dtype("int")), - ("column", np.dtype("int")), - ] + [ - (enforce_10ch_limit([name])[0], array_dict[name].dtype) - for name in names[3:] - ] - node = list(range(1, mg.ncol * mg.nrow + 1)) - col = list(range(1, mg.ncol + 1)) * mg.nrow - row = sorted(list(range(1, mg.nrow + 1)) * mg.ncol) - at = np.vstack( - [node, row, col] + [array_dict[name].ravel() for name in names[3:]] - ).transpose() - - names = enforce_10ch_limit(names) - - elif mg.grid_type == "vertex": - names = ["node"] + list(array_dict.keys()) - dtypes = [("node", np.dtype("int"))] + [ - (enforce_10ch_limit([name])[0], array_dict[name].dtype) - for name in names[1:] - ] - node = list(range(1, mg.ncpl + 1)) - at = np.vstack( - [node] + [array_dict[name].ravel() for name in names[1:]] - ).transpose() - - names = enforce_10ch_limit(names) - - elif mg.grid_type == "unstructured": - if mg.nlay is None: - names = ["node"] + list(array_dict.keys()) - dtypes = [("node", np.dtype("int"))] + [ - (enforce_10ch_limit([name])[0], array_dict[name].dtype) - for name in names[1:] - ] - node = list(range(1, mg.nnodes + 1)) - at = np.vstack( - [node] + [array_dict[name].ravel() for name in names[1:]] - ).transpose() + gdf = modelgrid.to_geodataframe() + names = list(array_dict.keys()) + at = np.vstack([array_dict[name].ravel() for name in names]) + names = enforce_10ch_limit(names) + + for ix, name in enumerate(names): + gdf[name] = at[ix] + + gdf = gdf.fillna(nan_val) + + if "epsg" in kwargs: + epsg = kwargs.pop("epsg") + crs = f"EPSG:{epsg}" + + if crs is not None: + if gdf.crs in None: + gdf = gdf.set_crs(crs) else: - names = ["node", "layer"] + list(array_dict.keys()) - dtypes = [("node", np.dtype("int")), ("layer", np.dtype("int"))] + [ - (enforce_10ch_limit([name])[0], array_dict[name].dtype) - for name in names[2:] - ] - node = list(range(1, mg.nnodes + 1)) - layer = np.zeros(mg.nnodes) - for ilay in range(mg.nlay): - istart, istop = mg.get_layer_node_range(ilay) - layer[istart:istop] = ilay + 1 - at = np.vstack( - [node] + [layer] + [array_dict[name].ravel() for name in names[2:]] - ).transpose() - - names = enforce_10ch_limit(names) - - # flag nan values and explicitly set the dtypes - if at.dtype in [float, np.float32, np.float64]: - at[np.isnan(at)] = nan_val - at = np.array([tuple(i) for i in at], dtype=dtypes) - - # write field information - fieldinfo = {name: get_pyshp_field_info(dtype.name) for name, dtype in dtypes} - for n in names: - w.field(n, *fieldinfo[n]) - - for i, r in enumerate(at): - # check if polygon is closed, if not close polygon for QGIS - if verts[i][-1] != verts[i][0]: - verts[i] = verts[i] + [verts[i][0]] - w.poly([verts[i]]) - w.record(*r) - - # close - w.close() + gdf = gdf.to_crs(crs) + + gdf.to_file(filename) + if verbose: - print(f"wrote {flopy_io.relpath_safe(path)}") + print(f"wrote {flopy_io.relpath_safe(os.getcwd(), filename)}") - # handle deprecated projection kwargs; warnings are raised in crs.py - write_prj_args = {} - if "epsg" in kwargs: - write_prj_args["epsg"] = kwargs.pop("epsg") - if "prj" in kwargs: - write_prj_args["prj"] = kwargs.pop("prj") - if kwargs: - raise TypeError(f"unhandled keywords: {kwargs}") - # write the projection file - try: - write_prj(path, mg, crs=crs, prjfile=prjfile, **write_prj_args) - except ImportError: - if verbose: - print("projection file not written") - return + if "prj" in kwargs or "prjfile" in kwargs or "wkt_string" in kwargs: + try: + write_prj(filename, modelgrid, crs=crs, prjfile=prjfile, **write_prj_args) + except ImportError: + if verbose: + print("projection file not written") def model_attributes_to_shapefile( - path: Union[str, PathLike], + filename: Union[str, PathLike], ml, package_names=None, array_dict=None, @@ -239,7 +161,7 @@ def model_attributes_to_shapefile( Parameters ---------- - path : str or PathLike + filename : str or PathLike path to write the shapefile to ml : flopy.mbase model instance @@ -276,6 +198,11 @@ def model_attributes_to_shapefile( >>> flopy.utils.model_attributes_to_shapefile('model.shp', m) """ + warnings.warn( + "model_attributes_to_shapefile is deprecated, please use the built in " + "to_geodataframe() method on the model object", + DeprecationWarning, + ) if array_dict is None: array_dict = {} @@ -286,143 +213,40 @@ def model_attributes_to_shapefile( else: package_names = [pak.name[0] for pak in ml.packagelist] - if "modelgrid" in kwargs: - grid = kwargs.pop("modelgrid") - else: - grid = ml.modelgrid - - horz_shape = grid.get_plottable_layer_shape() - for pname in package_names: - pak = ml.get_package(pname) - attrs = dir(pak) - if pak is not None: - if "start_datetime" in attrs: - attrs.remove("start_datetime") - for attr in attrs: - a = getattr(pak, attr) - if a is None or not hasattr(a, "data_type") or a.name == "thickness": - continue - if a.data_type == DataType.array2d: - if a.array is None or a.array.shape != horz_shape: - warn( - "Failed to get data for " - f"{a.name} array, {pak.name[0]} package" - ) - continue - name = shape_attr_name(a.name, keep_layer=True) - array_dict[name] = a.array - elif a.data_type == DataType.array3d: - # Not sure how best to check if an object has array data - if a.array is None: - warn( - "Failed to get data for " - f"{a.name} array, {pak.name[0]} package" - ) - continue - if isinstance(a.name, list) and a.name[0] == "thickness": - continue - if a.array.shape == horz_shape: - if hasattr(a, "shape"): - if a.shape[1] is None: # usg unstructured Util3d - # return a flattened array, - # with a.name[0] (a per-layer list) - array_dict[a.name[0]] = a.array - else: - array_dict[a.name] = a.array - else: - array_dict[a.name] = a.array - else: - # array is not the same shape as the layer shape - for ilay in range(a.array.shape[0]): - try: - arr = a.array[ilay] - except: - arr = a[ilay] - - if isinstance(a, Util3d): - aname = shape_attr_name(a[ilay].name) - else: - aname = a.name - - if arr.shape == (1,) + horz_shape: - # fix for mf6 case - arr = arr[0] - assert arr.shape == horz_shape - name = f"{aname}_{ilay + 1}" - array_dict[name] = arr - elif a.data_type == DataType.transient2d: - # Not sure how best to check if an object has array data - try: - assert a.array is not None - except: - warn( - "Failed to get data for " - f"{a.name} array, {pak.name[0]} package" - ) - continue - for kper in range(a.array.shape[0]): - name = f"{shape_attr_name(a.name)}{kper + 1}" - arr = a.array[kper][0] - assert arr.shape == horz_shape - array_dict[name] = arr - elif a.data_type == DataType.transientlist: - # Skip empty transientlist - if not a.data: - continue - - # Use first recarray kper to check transientlist - for kper in a.data.keys(): - if isinstance(a.data[kper], np.recarray): - break - # Skip transientlist if all elements are of object type - if all( - dtype == np.object_ - for dtype, _ in a.data[kper].dtype.fields.values() - ): - continue - - for name, array in a.masked_4D_arrays_itr(): - n = shape_attr_name(name, length=4) - for kper in range(array.shape[0]): - # guard clause for disu case - # array is (kper, node) only - if len(array.shape) == 2: - aname = f"{n}{kper + 1}" - arr = array[kper] - assert arr.shape == horz_shape - if np.all(np.isnan(arr)): - continue - array_dict[aname] = arr - continue - # non-disu case - for k in range(array.shape[1]): - aname = f"{n}{k + 1}{kper + 1}" - arr = array[kper][k] - assert arr.shape == horz_shape - if np.all(np.isnan(arr)): - continue - array_dict[aname] = arr - elif isinstance(a, list): - for v in a: - if ( - isinstance(a, DataInterface) - and v.data_type == DataType.array3d - ): - for ilay in range(a.model.modelgrid.nlay): - u2d = a[ilay] - name = f"{shape_attr_name(u2d.name)}_{ilay + 1}" - arr = u2d.array - assert arr.shape == horz_shape - array_dict[name] = arr - - # write data arrays to a shapefile - write_grid_shapefile(path, grid, array_dict) + gdf = ml.to_geodataframe(package_names=package_names, truncate_attrs=True) + + if array_dict: + modelgrid = ml.modelgrid + for name, array in array_dict.items(): + if modelgrid.grid_type() == "unstructured": + gdf[name] = array.ravel() + else: + if array.size == modelgrid.ncpl: + gdf[name] = array.ravel() + elif array.size % modelgrid.ncpl == 0: + array = array.reshape((-1, modelgrid.ncpl)) + for ix, lay in enumerate(array): + gdf[f"{name}_{ix}"] = lay + else: + raise ValueError( + f"{name} array size is not a multiple of ncpl {modelgrid.ncpl}" + ) + crs = kwargs.get("crs", None) + if crs is not None: + if gdf.crs is None: + gdf = gdf.set_crs(crs) + else: + gdf = gdf.to_crs(crs) + + gdf.to_file(filename) + prjfile = kwargs.get("prjfile", None) - try: - write_prj(path, grid, crs=crs, prjfile=prjfile) - except ImportError: - pass + if prjfile is not None: + try: + write_prj(filename, ml.modelgrid, crs=crs, prjfile=prjfile) + except ImportError: + pass def shape_attr_name(name, length=6, keep_layer=False): @@ -506,42 +330,12 @@ def truncate(s): return names -def get_pyshp_field_info(dtypename): - """Get pyshp dtype information for a given numpy dtype.""" - fields = { - "int": ("N", 18, 0), - "`, + such as an authority string (eg "EPSG:26916") or a WKT string. + **kwargs : Returns ------- @@ -1640,7 +1668,9 @@ def export_contours( from matplotlib.path import Path from ..utils.geometry import LineString - from .shapefile_utils import recarray2shp + from ..utils.utl_import import import_optional_dependency + + gpd = import_optional_dependency("geopandas") if not isinstance(contours, list): contours = [contours] @@ -1702,12 +1732,17 @@ def export_contours( if verbose: print(f"Writing {len(level)} contour lines") - ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray) + gdf = gpd.GeoDataFrame.from_features(geoms) + gdf[fieldname] = level + if crs is not None: + gdf = gdf.set_crs(crs) + gdf.to_file(filename) + return gdf - recarray2shp(ra, geoms, filename, **kwargs) - -def export_contourf(filename, contours, fieldname="level", verbose=False, **kwargs): +def export_contourf( + filename, contours, fieldname="level", verbose=False, crs=None, **kwargs +): """ Write matplotlib filled contours to shapefile. @@ -1723,11 +1758,17 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar the range represented by the polygon. Default is 'level'. verbose : bool, optional, default False whether to show verbose output + crs : pyproj.CRS, int, str, optional if `prjfile` is specified + Coordinate reference system (CRS) for the model grid + (must be projected; geographic CRS are not supported). + The value can be anything accepted by + :meth:`pyproj.CRS.from_user_input() `, + such as an authority string (eg "EPSG:26916") or a WKT string. **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp Returns ------- - None + gdf : GeoDataFrame of matplotlib contours Examples -------- @@ -1744,7 +1785,9 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar from matplotlib.path import Path from ..utils.geometry import Polygon, is_clockwise - from .shapefile_utils import recarray2shp + from ..utils.utl_import import import_optional_dependency + + gpd = import_optional_dependency("geopandas") if not isinstance(contours, list): contours = [contours] @@ -1829,9 +1872,12 @@ def export_contourf(filename, contours, fieldname="level", verbose=False, **kwar if verbose: print(f"Writing {len(level)} polygons") - ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray) - - recarray2shp(ra, geoms, filename, **kwargs) + gdf = gpd.GeoDataFrame.from_features(geoms) + gdf[fieldname] = level + if crs is not None: + gdf = gdf.set_crs(crs) + gdf.to_file(filename) + return gdf def export_array_contours( @@ -1842,6 +1888,7 @@ def export_array_contours( interval=None, levels=None, maxlevels=1000, + crs=None, **kwargs, ): """ @@ -1863,7 +1910,13 @@ def export_array_contours( list of contour levels maxlevels : int maximum number of contour levels - **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp + crs : pyproj.CRS, int, str, optional if `prjfile` is specified + Coordinate reference system (CRS) for the model grid + (must be projected; geographic CRS are not supported). + The value can be anything accepted by + :meth:`pyproj.CRS.from_user_input() `, + such as an authority string (eg "EPSG:26916") or a WKT string. + **kwargs : """ import matplotlib.pyplot as plt @@ -1880,7 +1933,7 @@ def export_array_contours( ctr = contour_array(modelgrid, ax, a, layer, levels=levels) kwargs["mg"] = modelgrid - export_contours(filename, ctr, fieldname, **kwargs) + export_contours(filename, ctr, fieldname, crs=crs, **kwargs) plt.close() diff --git a/flopy/mbase.py b/flopy/mbase.py index ef5a26b4be..f0c28f2b5c 100644 --- a/flopy/mbase.py +++ b/flopy/mbase.py @@ -760,6 +760,58 @@ def _output_msg(self, i, add=True): f"{txt2} the output list." ) + def to_geodataframe( + self, gdf=None, kper=0, package_names=None, truncate_attrs=False + ): + """ + Method to build a Geodataframe from model inputs. Note: transient data + will only be exported for a single stress period. + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + package_names : list + optional list of package names + truncate_attrs : bool + method to truncate attribute names for shapefile attribute name length + restrictions + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + modelgrid = self.modelgrid + if modelgrid is not None: + gdf = modelgrid.to_geodataframe() + else: + raise AttributeError( + "model does not have a grid instance, please supply a geodataframe" + ) + + # todo: needs testing + print("break") + if package_names is None: + package_names = [pak.name[0] for pak in self.packagelist] + else: + pass + + for package in self.packagelist: + if package.package_type in ("hfb6",): + continue + if package.name[0] not in package_names: + # todo: this needs testing + continue + if callable(getattr(package, "to_geodataframe", None)): + gdf = package.to_geodataframe( + gdf, kper=kper, sparse=False, truncate_attrs=truncate_attrs + ) + + return gdf + def add_output_file( self, unit, diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py index e00c877b1a..dc9b70ea3f 100644 --- a/flopy/mf6/data/mfdataarray.py +++ b/flopy/mf6/data/mfdataarray.py @@ -323,6 +323,71 @@ def data(self): """Returns array data. Calls get_data with default parameters.""" return self._get_data() + def to_geodataframe(self, gdf=None, name=None, forgive=False, truncate_attrs=False, **kwargs): + """ + Method to add an input array to a geopandas GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame object + name : str + optional attribute name, default uses util2d name + forgive : bool + optional flag to continue if data shape not compatible with GeoDataFrame + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + geopandas GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.to_geodataframe() + + if modelgrid is not None: + if modelgrid.grid_type != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + if name is None: + name = self.name + + data = self.array + if data is None: + return gdf + + if truncate_attrs: + name = shape_attr_name(name=name) + + if data.size == ncpl: + gdf[name] = data.ravel() + + elif data.size % ncpl == 0: + data = data.reshape((-1, ncpl)) + for ix, arr in enumerate(data): + aname = f"{name}_{ix}" + gdf[aname] = arr + elif forgive: + return gdf + else: + raise ValueError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def new_simulation(self, sim_data): """Initialize MFArray object for a new simulation @@ -1890,6 +1955,76 @@ def _build_period_data( output[sp] = data return output + def to_geodataframe(self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs): + """ + Method to add an input array to a geopandas GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame object + name : str + optional attribute name, default uses util2d name + kper : int + stress period number + forgive : bool + optional flag to continue if data shape not compatible with GeoDataFrame + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + geopandas GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.to_geodataframe() + + if modelgrid is not None: + if modelgrid.grid_type != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + if self.array is None: + return gdf + + if truncate_attrs: + name = shape_attr_name(self.name, length=4) + else: + name = f"{self.path[1]}_{self.name}" + + data = self.get_data(key=kper, apply_mult=True) + if data.size == ncpl: + name = f"{name}_{kper}" + gdf[name] = data.ravel() + + elif data.size % ncpl == 0: + data = data.reshape((-1, ncpl)) + for ix, arr in enumerate(data): + if truncate_attrs: + aname = f"{name}{ix}{kper}" + else: + aname = f"{name}_{ix}_{kper}" + gdf[aname] = arr + elif forgive: + return gdf + else: + raise ValueError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def set_record(self, data_record): """Sets data and metadata at layer `layer` and time `key` to `data_record`. For unlayered data do not pass in `layer`. diff --git a/flopy/mf6/data/mfdatalist.py b/flopy/mf6/data/mfdatalist.py index 529177c41c..d089caa747 100644 --- a/flopy/mf6/data/mfdatalist.py +++ b/flopy/mf6/data/mfdatalist.py @@ -145,6 +145,63 @@ def to_array(self, kper=0, mask=False): model_grid = self.data_dimensions.get_model_grid() return list_to_array(sarr, model_grid, kper, mask) + def to_geodataframe(self, gdf=None, sparse=False, truncate_attrs=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + data = self.to_array(mask=True) + if data is None: + return gdf + + col_names = [] + for name, array3d in data.items(): + if truncate_attrs: + aname = shape_attr_name(name) + else: + aname = f"{self.path[1].lower()}_{name}" + + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def new_simulation(self, sim_data): """Initialize MFList object for a new simulation. @@ -1596,6 +1653,69 @@ def to_array(self, kper=0, mask=False): """Returns list data as an array.""" return super().to_array(kper, mask) + def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + data = self.to_array(kper=kper, mask=True) + if data is None: + return gdf + + col_names = [] + for name, array3d in data.items(): + if truncate_attrs: + name = shape_attr_name(name, length=4) + else: + name = f"{self.path[1].lower()}_{name}" + + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[name] = array + col_names.append(name) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + if truncate_attrs: + aname = f"{name}{lay}{kper}" + else: + aname = f"{name}_{lay}_{kper}" + gdf[aname] = arr.ravel() + col_names.append(aname) + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def remove_transient_key(self, transient_key): """Remove transient stress period key. Method is used internally by FloPy and is not intended to the end user. diff --git a/flopy/mf6/data/mfdataplist.py b/flopy/mf6/data/mfdataplist.py index a7aef4625e..bb6c1e2af0 100644 --- a/flopy/mf6/data/mfdataplist.py +++ b/flopy/mf6/data/mfdataplist.py @@ -842,6 +842,63 @@ def to_array(self, kper=0, mask=False): model_grid = self.data_dimensions.get_model_grid() return list_to_array(sarr, model_grid, kper, mask) + def to_geodataframe(self, gdf=None, sparse=False, truncate_attrs=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + data = self.to_array(mask=True) + if data is None: + return gdf + + col_names = [] + for name, array3d in data.items(): + if truncate_attrs: + aname = shape_attr_name(name) + else: + aname = f"{self.path[1].lower()}_{name}" + + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + gdf[f"{aname}_{lay}"] = arr.ravel() + col_names.append(f"{aname}_{lay}") + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def set_record(self, record, autofill=False, check_data=True): """Sets the contents of the data and metadata to "data_record". Data_record is a dictionary with has the following format: @@ -2019,6 +2076,68 @@ def __init__( self.repeating = True self.empty_keys = {} + def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + from ...export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + data = self.to_array(kper=kper, mask=True) + + col_names = [] + for name, array3d in data.items(): + if truncate_attrs: + name = shape_attr_name(name, length=4) + else: + name = f"{self.path[1].lower()}_{name}" + if modelgrid.grid_type == "unstructured": + array = array3d.ravel() + aname = f"{name}_{kper}" + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array3d[lay].ravel() + if truncate_attrs: + aname = f"{name}{lay}{kper}" + else: + aname = f"{name}_{lay}_{kper}" + gdf[aname] = arr.ravel() + col_names.append(aname) + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + + @property def data_type(self): return DataType.transientlist diff --git a/flopy/mf6/data/mfdatautil.py b/flopy/mf6/data/mfdatautil.py index 9b0a885119..17ee03ab8c 100644 --- a/flopy/mf6/data/mfdatautil.py +++ b/flopy/mf6/data/mfdatautil.py @@ -202,6 +202,8 @@ def list_to_array(sarr, model_grid, kper=0, mask=False): cnt = np.zeros(shape, dtype=np.float64) for sp_rec in sarr: if sp_rec is not None: + if "cellid" not in sp_rec.dtype.names: + return None for rec in sp_rec: arr[rec["cellid"]] += rec[name] cnt[rec["cellid"]] += 1.0 diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py index ab2a845f44..2282ab1583 100644 --- a/flopy/mf6/mfmodel.py +++ b/flopy/mf6/mfmodel.py @@ -793,6 +793,57 @@ def output(self): except AttributeError: return MF6Output(self, budgetkey=budgetkey) + def to_geodataframe(self, gdf=None, kper=0, package_names=None, truncate_attrs=False): + """ + Method to build a Geodataframe from model inputs. Note: transient data + will only be exported for a single stress period. + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + package_names : list + optional list of package names + truncate_attrs : bool + method to truncate attribute names for shapefile attribute name length + restrictions + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + modelgrid = self.modelgrid + if modelgrid is not None: + gdf = modelgrid.to_geodataframe() + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + + print('break') + # todo: need to check this one + if package_names is None: + package_names = [pak.name[0] for pak in self.packagelist] + else: + pass + + for package in self.packagelist: + if package.package_type in ("hfb",): + continue + if package.name not in package_names and package.package_type not in package_names: + # todo: this needs testing + continue + if callable(getattr(package, "to_geodataframe", None)): + gdf = package.to_geodataframe( + gdf, kper=kper, sparse=False, truncate_attrs=truncate_attrs + ) + + return gdf + def export(self, f, **kwargs): """Method to export a model to a shapefile or netcdf file diff --git a/flopy/mf6/mfpackage.py b/flopy/mf6/mfpackage.py index 702becfccd..98166e5734 100644 --- a/flopy/mf6/mfpackage.py +++ b/flopy/mf6/mfpackage.py @@ -13,6 +13,7 @@ from ..pakbase import PackageInterface from ..utils import datautil from ..utils.check import mf6check +from ..utils.utl_import import import_optional_dependency from ..version import __version__ from .coordinates import modeldimensions from .data import ( @@ -2088,6 +2089,84 @@ def package_filename_dict(self): ) return self._package_container.package_filename_dict + def to_geodataframe(self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs): + """ + Method to create a GeoDataFrame from a modflow package + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + sparse : bool + optional parameter to create a sparse geodataframe + truncate_attrs : bool + method to truncate attribute names for shapefile attribute name length + restrictions + + Returns + ------- + gdf : GeoDataFrame + """ + if gdf is None: + if isinstance(self.parent, ModelInterface): + modelgrid = self.parent.modelgrid + if modelgrid is not None: + if self.package_type == "hfb": + gpd = import_optional_dependency("geopandas") + from ..plot.plotutil import hfb_data_to_linework + + recarray = self.stress_period_data.data[kper] + lines = hfb_data_to_linework(recarray, modelgrid) + geo_interface = {"type": "FeatureCollection"} + features = [ + { + "id": f"{ix}", + "geometry": {"coordinates": line, "type": "LineString"}, + "properties": {} + } + for ix, line in enumerate(lines) + ] + geo_interface["features"] = features + gdf = gpd.GeoDataFrame.from_features(geo_interface) + + for name in recarray.dtype.names: + gdf[name] = recarray[name] + + return gdf + + else: + gdf = modelgrid.to_geodataframe() + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + else: + raise AssertionError( + "Package does not have a model instance, " + "please supply a geodataframe" + ) + + for attr, value in self.__dict__.items(): + if callable(getattr(value, "to_geodataframe", None)): + if isinstance(value, (ModelInterface, PackageInterface)): + continue + # do not pass sparse in here, "sparsify" after all data has been + # added to geodataframe + gdf = value.to_geodataframe( + gdf, forgive=True, kper=kper, sparse=False, truncate_attrs=truncate_attrs + ) + + if sparse: + col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + + def get_package(self, name=None, type_only=False, name_only=False): """ Finds a package by package name, package key, package type, or partial diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py index 468cf1a421..70cb7cae4a 100644 --- a/flopy/modflow/mfhfb.py +++ b/flopy/modflow/mfhfb.py @@ -16,6 +16,7 @@ from ..pakbase import Package from ..utils.flopy_io import line_parse from ..utils.recarray_utils import create_empty_recarray +from ..utils.utl_import import import_optional_dependency from .mfparbc import ModflowParBc as mfparbc @@ -181,6 +182,47 @@ def _ncells(self): """ return self.nhfbnp + def _get_hfb_lines(self): + """ + Method to get hfb lines. Lines are ordered by recarray index + + Returns + ------- + list : list of hfb lines + """ + from ..plot.plotutil import hfb_data_to_linework + + return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid) + + def to_geodataframe(self, **kwargs): + """ + Method to create a LineString based geodataframe of horizontal flow barriers + + Returns + ------- + GeoDataFrame + + """ + gpd = import_optional_dependency("geopandas") + + lines = self._get_hfb_lines() + geo_interface = {"type": "FeatureCollection"} + features = [ + { + "id": f"{ix}", + "geometry": {"coordinates": line, "type": "LineString"}, + "properties": {}, + } + for ix, line in enumerate(lines) + ] + geo_interface["features"] = features + gdf = gpd.GeoDataFrame.from_features(geo_interface) + hfbs = self.hfb_data + for name in hfbs.dtype.names: + gdf[name] = hfbs[name] + + return gdf + def write_file(self): """ Write the package file. diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index 306d77a661..c34f59d7d1 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -622,6 +622,43 @@ def paths(self): def df(self): return pd.DataFrame(self.reach_data) + def to_geodataframe(self, gdf=None, sparse=True, **kwargs): + """ + Method to export SFR reach data to a GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + sparse : bool + boolean flag for sparse dataframe construction. Default is True + """ + modelgrid = self.parent.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + df = self.df + if "k" in list(df): + df["node"] = df["node"] - (modelgrid.ncpl * df["k"]) + df = df.drop(columns=["k", "i", "j"]) + + df["node"] += 1 + gdf = gdf.merge(df, how="left", on="node") + + if sparse: + col_names = [col for col in list(df) if col != "node"] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + else: + gdf = gdf.drop_duplicates(subset=["node"]) + gdf = gdf.reset_index(drop=True) + + return gdf + def _make_graph(self): # get all segments and their outseg graph = {} diff --git a/flopy/pakbase.py b/flopy/pakbase.py index 70c6a23ddc..0a2c5ed26c 100644 --- a/flopy/pakbase.py +++ b/flopy/pakbase.py @@ -674,6 +674,67 @@ def export(self, f, **kwargs): return export.utils.package_export(f, self, **kwargs) + def to_geodataframe( + self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs + ): + """ + Method to create a GeoDataFrame from a modflow package + + Parameters + ---------- + gdf : GeoDataFrame + optional geopandas geodataframe object to add data to. Default is None + kper : int + stress period to get transient data from + sparse : bool + optional parameter to create a sparse geodataframe + truncate_attrs : bool + method to truncate attribute names for shapefile attribute name length + restrictions + + Returns + ------- + gdf : GeoDataFrame + """ + from .mbase import BaseModel + + if gdf is None: + if isinstance(self.parent, BaseModel): + modelgrid = self.parent.modelgrid + if modelgrid is not None: + gdf = modelgrid.to_geodataframe() + else: + raise AttributeError( + "model does not have a grid instance, " + "please supply a geodataframe" + ) + else: + raise AssertionError( + "Package does not have a model instance, " + "please supply a geodataframe" + ) + + for attr, value in self.__dict__.items(): + if callable(getattr(value, "to_geodataframe", None)): + if isinstance(value, (BaseModel, PackageInterface)): + continue + # do not pass sparse in here, make sparse after all data has been + # added to geodataframe + gdf = value.to_geodataframe( + gdf, + forgive=True, + kper=kper, + sparse=False, + truncate_attrs=truncate_attrs, + ) + + if sparse: + col_names = [i for i in gdf if i not in ("geometry", "node", "row", "col")] + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def _generate_heading(self): """Generate heading.""" from . import __version__ diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py index e8085cc8f2..d6ae28f694 100644 --- a/flopy/plot/plotutil.py +++ b/flopy/plot/plotutil.py @@ -2947,3 +2947,71 @@ def to_prt_pathlines( return df else: return ret + + +def hfb_data_to_linework(recarray, modelgrid): + """ + Method to get lines representing horizontal flow barriers + + Parameters + ---------- + recarray : np.recarray + recarray of hfb input data + modelgrid : Grid + modelgrid instance + + Returns + ------- + list : list of line segments + """ + iverts = modelgrid.iverts + verts = modelgrid.verts + nodes = [] + if modelgrid.grid_type == "structured": + if "k" in recarray.dtype.names: + for rec in recarray: + node1 = modelgrid.get_node([(0, rec["irow1"], rec["icol1"])])[0] + node2 = modelgrid.get_node([(0, rec["irow2"], rec["icol2"])])[0] + nodes.append((node1, node2)) + else: + for rec in recarray: + node1 = modelgrid.get_node([(0,) + rec["cellid1"][1:]])[0] + node2 = modelgrid.get_node([(0,) + rec["cellid2"][1:]])[0] + nodes.append((node1, node2)) + + elif modelgrid.grid_type == "vertex": + for rec in recarray: + nodes.append((rec["cellid1"][-1], rec["cellid2"][-1])) + + else: + if "node1" in recarray.dtype.names: + nodes = list(zip(recarray["node1"], recarray["node2"])) + else: + for rec in recarray: + nodes.append((rec["cellid1"][0], rec["cellid2"][0])) + + shared_edges = [] + for node0, node1 in nodes: + iv0 = iverts[node0] + iv1 = iverts[node1] + edges = [] + for ix in range(len(iv0)): + edges.append(tuple(sorted((iv0[ix - 1], iv0[ix])))) + + for ix in range(len(iv1)): + edge = tuple(sorted((iv1[ix - 1], iv1[ix]))) + if edge in edges: + shared_edges.append(edge) + break + + if not shared_edges: + raise AssertionError( + f"No shared cell edges found. Cannot represent HFB " + f"for nodes {node0} and {node1}" + ) + + lines = [] + for edge in shared_edges: + lines.append((tuple(verts[edge[0]]), tuple(verts[edge[1]]))) + + return lines diff --git a/flopy/utils/binaryfile/__init__.py b/flopy/utils/binaryfile/__init__.py index f0781d56a3..2cfd79a863 100644 --- a/flopy/utils/binaryfile/__init__.py +++ b/flopy/utils/binaryfile/__init__.py @@ -1604,6 +1604,56 @@ def get_data( if t ] + def to_geodataframe( + self, + gdf=None, + modelgrid=None, + idx=None, + kstpkper=None, + totim=None, + text=None, + ): + if (idx, kstpkper, totim) == (None, None, None): + raise AssertionError( + "to_geodataframe() missing 1 required argument: " + "please provide 'idx', 'kstpkper', or 'totim'" + ) + + if gdf is None: + if modelgrid is None: + if self.modelgrid is None: + raise AssertionError( + "A geodataframe or modelgrid instance must be supplied" + ) + modelgrid = self.modelgrid + + gdf = modelgrid.to_geodataframe() + + col_names = [] + if text is None: + textlist = [i.decode() for i in self.textlist] + for text in textlist: + data = self.get_data( + idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True + ) + + text = text.strip().lower().replace(" ", "_") + for ix, arr in enumerate(data[0]): + name = f"{text}_{ix}" + gdf[name] = arr.ravel() + col_names.append(name) + else: + data = self.get_data( + idx=idx, kstpkper=kstpkper, totim=totim, text=text, full3D=True + ) + text = text.strip().lower().replace(" ", "_") + for ix, arr in enumerate(data[0]): + name = f"{text}_{ix}" + gdf[name] = arr.ravel() + col_names.append(name) + + return gdf + def get_ts(self, idx, text=None, times=None, variable="q"): """ Get a time series from the binary budget file. diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index 6f81268fed..845cc73e5b 100644 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -197,18 +197,18 @@ def __init__(self, filename: Union[str, PathLike], precision, verbose, **kwargs) self.model = None self.dis = None - self.mg = None + self.modelgrid = None if "model" in kwargs.keys(): self.model = kwargs.pop("model") - self.mg = self.model.modelgrid + self.modelgrid = self.model.modelgrid self.dis = self.model.dis if "dis" in kwargs.keys(): self.dis = kwargs.pop("dis") - self.mg = self.dis.parent.modelgrid + self.modelgrid = self.dis.parent.modelgrid if "tdis" in kwargs.keys(): self.tdis = kwargs.pop("tdis") if "modelgrid" in kwargs.keys(): - self.mg = kwargs.pop("modelgrid") + self.modelgrid = kwargs.pop("modelgrid") if len(kwargs.keys()) > 0: args = ",".join(kwargs.keys()) raise Exception(f"LayerFile error: unrecognized kwargs: {args}") @@ -217,9 +217,9 @@ def __init__(self, filename: Union[str, PathLike], precision, verbose, **kwargs) self._build_index() # now that we read the data and know nrow and ncol, - # we can make a generic mg if needed - if self.mg is None: - self.mg = StructuredGrid( + # we can make a generic modelgrid if needed + if self.modelgrid is None: + self.modelgrid = StructuredGrid( delc=np.ones((self.nrow,)), delr=np.ones(self.ncol), nlay=self.nlay, @@ -240,6 +240,54 @@ def __enter__(self): def __exit__(self, *exc): self.close() + def to_geodataframe( + self, gdf=None, modelgrid=None, kstpkper=None, totim=None, attrib_name=None + ): + """ + Generate a GeoDataFrame with data from a LayerFile instance + + Parameters + ---------- + gdf : GeoDataFrame + optional, existing geodataframe with NCPL geometries + modelgrid : Grid + optional modelgrid instance to generate a GeoDataFrame from + kstpkper : tuple of ints + A tuple containing the time step and stress period (kstp, kper). + These are zero-based kstp and kper values. + totim : float + The simulation time. + attrib_name : str + optional base name of attribute columns. (default is text attribute) + + + Returns + ------- + GeoDataFrame + """ + if gdf is None: + if modelgrid is None: + if self.modelgrid is None: + raise AssertionError( + "A geodataframe or modelgrid instance must be supplied" + ) + modelgrid = self.modelgrid + + gdf = modelgrid.to_geodataframe() + + array = np.atleast_3d( + self.get_data(kstpkper=kstpkper, totim=totim).transpose() + ).transpose() + + if attrib_name is None: + attrib_name = self.text.decode() + + for ix, arr in enumerate(array): + name = f"{attrib_name}_{ix}" + gdf[name] = np.ravel(arr) + + return gdf + def to_shapefile( self, filename: Union[str, PathLike], @@ -287,7 +335,10 @@ def to_shapefile( >>> times = hdobj.get_times() >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1]) """ - + warnings.warn( + "to_shapefile() is deprecated and is being replaced by to_geodataframe()", + DeprecationWarning, + ) plotarray = np.atleast_3d( self.get_data(kstpkper=kstpkper, totim=totim, mflay=mflay).transpose() ).transpose() @@ -301,7 +352,7 @@ def to_shapefile( from ..export.shapefile_utils import write_grid_shapefile - write_grid_shapefile(filename, self.mg, attrib_dict, verbose) + write_grid_shapefile(filename, self.modelgrid, attrib_dict, verbose) def plot( self, @@ -413,7 +464,7 @@ def plot( axes=axes, filenames=filenames, mflay=mflay, - modelgrid=self.mg, + modelgrid=self.modelgrid, **kwargs, ) diff --git a/flopy/utils/geospatial_utils.py b/flopy/utils/geospatial_utils.py index a60a16aca5..36cfde7b1e 100644 --- a/flopy/utils/geospatial_utils.py +++ b/flopy/utils/geospatial_utils.py @@ -453,7 +453,7 @@ def shapely(self): return self._shapely @property - def geo_dataframe(self): + def geodataframe(self): """ Property that returns a geopandas DataFrame diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index a5bc077886..5a34175d5f 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -12,6 +12,7 @@ from flopy.utils.particletrackfile import ParticleTrackFile from ..utils.flopy_io import loadtxt +from ..utils.utl_import import import_optional_dependency class ModpathFile(ParticleTrackFile): @@ -671,6 +672,61 @@ def get_destination_endpoint_data(self, dest_cells, source=False): inds = np.isin(raslice, dest_cells) return data[inds].copy().view(np.recarray) + def to_geodataframe( + self, + modelgrid, + data=None, + direction="ending", + ): + """ + Create a geodataframe of particle starting / ending locations. + + Parameters + ---------- + modelgrid : flopy.discretization.grid instance + Used to scale and rotate Global x,y,z values in MODPATH Endpoint + file. + data : np.recarray + Record array of same form as that returned by EndpointFile.get_alldata. + (if none, EndpointFile.get_alldata() is exported). + direction : str + String defining if 'starting' or 'ending' particle locations should be + considered. (default is 'ending') + """ + from ..utils import geometry + + gpd = import_optional_dependency("geopandas") + shapely_geo = import_optional_dependency("shapely.geometry") + if data is None: + data = self.get_alldata() + + if direction.lower() == "ending": + xcol, ycol, zcol = "x", "y", "z" + elif direction.lower() == "starting": + xcol, ycol, zcol = "x0", "y0", "z0" + else: + raise Exception( + 'flopy.map.plot_endpoint direction must be "ending" or "starting".' + ) + x, y = geometry.transform( + data[xcol], + data[ycol], + xoff=modelgrid.xoffset, + yoff=modelgrid.yoffset, + angrot_radians=modelgrid.angrot_radians, + ) + z = data[zcol] + + geoms = [shapely_geo.Point(p) for p in zip(x, y, z)] + gdf = gpd.GeoDataFrame(data, geometry=geoms, crs=modelgrid.crs) + + # adjust to 1 based node numbers + for col in list(gdf): + if col in self.kijnames: + gdf[col] += 1 + + return gdf + def write_shapefile( self, data=None, @@ -715,44 +771,21 @@ def write_shapefile( - ``epsg`` (int): use ``crs`` instead. """ - from ..discretization import StructuredGrid - from ..export.shapefile_utils import recarray2shp - from ..utils import geometry - from ..utils.geometry import Point + import warnings + warnings.warn( + "write_shapefile is Deprecated, please use to_geodataframe() in the future" + ) epd = (data if data is not None else endpoint_data).copy() - if epd is None: - epd = self.get_alldata() + gdf = self.to_geodataframe(modelgrid=mg, data=epd, direction=direction) - if direction.lower() == "ending": - xcol, ycol, zcol = "x", "y", "z" - elif direction.lower() == "starting": - xcol, ycol, zcol = "x0", "y0", "z0" - else: - raise Exception( - 'flopy.map.plot_endpoint direction must be "ending" or "starting".' - ) - if mg is None: - raise ValueError("A modelgrid object was not provided.") - - if isinstance(mg, StructuredGrid): - x, y = geometry.transform( - epd[xcol], - epd[ycol], - xoff=mg.xoffset, - yoff=mg.yoffset, - angrot_radians=mg.angrot_radians, - ) - else: - x, y = mg.get_coords(epd[xcol], epd[ycol]) - z = epd[zcol] + if crs is not None: + if gdf.crs is None: + gdf = gdf.set_crs(crs) + else: + gdf = gdf.to_crs(crs) - geoms = [Point(x[i], y[i], z[i]) for i in range(len(epd))] - # convert back to one-based - for n in self.kijnames: - if n in epd.dtype.names: - epd[n] += 1 - recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs) + gdf.to_file(shpname) class TimeseriesFile(ModpathFile): diff --git a/flopy/utils/particletrackfile.py b/flopy/utils/particletrackfile.py index 00e9a7a141..47db7f89ef 100644 --- a/flopy/utils/particletrackfile.py +++ b/flopy/utils/particletrackfile.py @@ -10,6 +10,8 @@ import numpy as np from numpy.lib.recfunctions import stack_arrays +from .utl_import import import_optional_dependency + MIN_PARTICLE_TRACK_DTYPE = np.dtype( [ ("x", np.float32), @@ -183,6 +185,125 @@ def intersect(self, cells, to_recarray) -> np.recarray: """Find intersection of pathlines with cells.""" pass + def to_geodataframe( + self, + modelgrid, + data=None, + one_per_particle=True, + direction="ending", + ): + """ + Create a geodataframe of particle tracks. + + Parameters + ---------- + modelgrid : flopy.discretization.Grid instance + Used to scale and rotate Global x, y, z values. + data : np.recarray + Record array of same form as that returned by + get_alldata(). (if none, get_alldata() is exported). + one_per_particle : boolean (default True) + True writes a single LineString with a single set of attribute + data for each particle. False writes a record/geometry for each + pathline segment (each row in the PathLine file). This option can + be used to visualize attribute information (time, model layer, + etc.) across a pathline in a GIS. + direction : str + String defining if starting or ending particle locations should be + included in shapefile attribute information. Only used if + one_per_particle=False. (default is 'ending') + + Returns + ------- + GeoDataFrame + """ + from . import geometry + + shapely_geo = import_optional_dependency("shapely.geometry") + gpd = import_optional_dependency("geopandas") + + if data is None: + data = self._data.view(np.recarray) + else: + # convert pathline list to a single recarray + if isinstance(data, list): + s = data[0] + print(s.dtype) + for n in range(1, len(data)): + s = stack_arrays((s, data[n])) + data = s.view(np.recarray) + + data = data.copy() + data.sort(order=["particleid", "time"]) + + particles = np.unique(data.particleid) + geoms = [] + + # create a dict of attrs? + headings = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] + attrs = [] + for h in headings: + if h in data.dtype.names: + attrs.append(h) + + if one_per_particle: + dfdata = {a: [] for a in attrs} + if direction == "ending": + idx = -1 + else: + idx = 0 + + for p in particles: + ra = data[data.particleid == p] + for k in dfdata.keys(): + if k == "time": + dfdata[k].append(np.max(ra[k])) + else: + dfdata[k].append(ra[k][idx]) + + x, y = geometry.transform( + ra.x, + ra.y, + modelgrid.xoffset, + modelgrid.yoffset, + modelgrid.angrot_radians, + ) + z = ra.z + + line = list(zip(x, y, z)) + geoms.append(shapely_geo.LineString(line)) + + else: + dfdata = {a: [] for a in data.dtype.names} + for pid in particles: + ra = data[data.particleid == pid] + x, y = geometry.transform( + ra.x, + ra.y, + modelgrid.xoffset, + modelgrid.yoffset, + modelgrid.angrot_radians, + ) + z = ra.z + geoms += [ + shapely_geo.LineString( + [(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])] + ) + for i in np.arange(1, (len(ra))) + ] + for k in dfdata.keys(): + dfdata[k].extend(ra[k][1:]) + + # now create a geodataframe + gdf = gpd.GeoDataFrame(dfdata, geometry=geoms, crs=modelgrid.crs) + + # adjust to 1 based node numbers + for col in list(gdf): + if col in self.kijnames: + gdf[col] += 1 + + return gdf + def write_shapefile( self, data=None, @@ -227,102 +348,22 @@ def write_shapefile( - ``epsg`` (int): use ``crs`` instead. """ - from ..export.shapefile_utils import recarray2shp - from . import geometry - from .geometry import LineString - - series = data - if series is None: - series = self._data.view(np.recarray) - else: - # convert pathline list to a single recarray - if isinstance(series, list): - s = series[0] - print(s.dtype) - for n in range(1, len(series)): - s = stack_arrays((s, series[n])) - series = s.view(np.recarray) - - series = series.copy() - series.sort(order=["particleid", "time"]) + import warnings - if mg is None: - raise ValueError("A modelgrid object was not provided.") - - particles = np.unique(series.particleid) - geoms = [] - - # create dtype with select attributes in pth - names = series.dtype.names - dtype = [] - atts = ["particleid", "particlegroup", "time", "k", "i", "j", "node"] - for att in atts: - if att in names: - t = np.int32 - if att == "time": - t = np.float32 - dtype.append((att, t)) - dtype = np.dtype(dtype) - - # reset names to the selected names in the created dtype - names = dtype.names - - # 1 geometry for each path - if one_per_particle: - loc_inds = 0 - if direction == "ending": - loc_inds = -1 - - sdata = [] - for pid in particles: - ra = series[series.particleid == pid] - - x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians - ) - z = ra.z - geoms.append(LineString(list(zip(x, y, z)))) - - t = [pid] - if "particlegroup" in names: - t.append(ra.particlegroup[0]) - t.append(ra.time.max()) - if "k" in names: - t.append(ra.k[loc_inds]) - if "node" in names: - t.append(ra.node[loc_inds]) - else: - if "i" in names: - t.append(ra.i[loc_inds]) - if "j" in names: - t.append(ra.j[loc_inds]) - sdata.append(tuple(t)) - - sdata = np.array(sdata, dtype=dtype).view(np.recarray) - - # geometry for each row in PathLine file - else: - dtype = series.dtype - sdata = [] - for pid in particles: - ra = series[series.particleid == pid] - if mg is not None: - x, y = geometry.transform( - ra.x, ra.y, mg.xoffset, mg.yoffset, mg.angrot_radians - ) - else: - x, y = geometry.transform(ra.x, ra.y, 0, 0, 0) - z = ra.z - geoms += [ - LineString([(x[i - 1], y[i - 1], z[i - 1]), (x[i], y[i], z[i])]) - for i in np.arange(1, (len(ra))) - ] - sdata += ra[1:].tolist() - sdata = np.array(sdata, dtype=dtype).view(np.recarray) - - # convert back to one-based - for n in set(self.kijnames).intersection(set(sdata.dtype.names)): - sdata[n] += 1 + warnings.warn( + "write_shapefile will be Deprecated, please use to_geo_dataframe()", + DeprecationWarning, + ) + gdf = self.to_geodataframe( + modelgrid=mg, + data=data, + one_per_particle=one_per_particle, + direction=direction, + ) + if crs is not None: + if gdf.crs is None: + gdf = gdf.set_crs(crs) + else: + gdf = gdf.to_crs(crs) - # write the final recarray to a shapefile - recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs) + gdf.to_file(shpname) diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py index 3823260615..afef0d8567 100644 --- a/flopy/utils/util_array.py +++ b/flopy/utils/util_array.py @@ -634,6 +634,55 @@ def data_type(self): def plottable(self): return True + def to_geodataframe( + self, gdf=None, forgive=False, name=None, truncate_attrs=False, **kwargs + ): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + forgive : bool + optional flag to continue running and pass data that is not compatible + with the geodataframe shape + name : str + optional array name base. If None, method uses the .name attribute + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + if name is not None: + names = [name for _ in range(len(self.util_2ds))] + else: + names = self.name + for lay, u2d in enumerate(self.util_2ds): + name = f"{names[lay]}_{lay}" + gdf = u2d.to_geodataframe( + gdf=gdf, + name=name, + forgive=forgive, + truncate_attrs=truncate_attrs, + **kwargs, + ) + + return gdf + def export(self, f, **kwargs): from .. import export @@ -1059,6 +1108,37 @@ def data_type(self): def plottable(self): return False + def to_geodataframe( + self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs + ): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + forgive : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + u3d = self.transient_3ds[kper] + # note: may need to provide a pass through name for u3d to avoid s.p. + # number being tacked on. Test this on a model with the LAK package... + name = self.name_base[:-1].lower() + gdf = u3d.to_geodataframe( + gdf=gdf, forgive=forgive, name=name, truncate_attrs=truncate_attrs, **kwargs + ) + return gdf + def get_zero_3d(self, kper): name = f"{self.name_base}{kper + 1}(filled zero)" return Util3d( @@ -1399,6 +1479,35 @@ def from_4d(cls, model, pak_name, m4ds): name=name, ) + def to_geodataframe( + self, gdf=None, kper=0, forgive=False, truncate_attrs=False, **kwargs + ): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + forgive : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + u2d = self.transient_2ds[kper] + name = self.name_base[:-1] + gdf = u2d.to_geodataframe( + gdf=gdf, name=name, forgive=forgive, truncate_attrs=truncate_attrs, **kwargs + ) + return gdf + def __setattr__(self, key, value): if hasattr(self, "transient_2ds") and key == "cnstnt": # set cnstnt for each u2d @@ -1876,6 +1985,65 @@ def _decide_how(self): else: self._how = "internal" + def to_geodataframe( + self, gdf=None, name=None, forgive=False, truncate_attrs=False, **kwargs + ): + """ + Method to add an input array to a geopandas GeoDataFrame + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame object + name : str + optional attribute name, default uses util2d name + forgive : bool + optional flag to continue if data shape not compatible with GeoDataFrame + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + geopandas GeoDataFrame + """ + from ..export.shapefile_utils import shape_attr_name + + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if gdf is None: + if modelgrid is None: + return gdf + gdf = modelgrid.to_geodataframe() + + if modelgrid is not None: + if modelgrid.grid_type != "unstructured": + ncpl = modelgrid.ncpl + else: + ncpl = modelgrid.nnodes + else: + ncpl = len(gdf) + + if name is None: + name = self.name + + if truncate_attrs: + name = shape_attr_name(name, keep_layer=True) + + data = self.array + + if data.size == ncpl: + gdf[name] = data.ravel() + elif forgive: + return gdf + else: + raise AssertionError( + f"Data size {data.size} not compatible with dataframe length {ncpl}" + ) + + return gdf + def plot( self, title=None, diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index 3ddbb550fa..c078c5cca5 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -14,6 +14,7 @@ import numpy as np import pandas as pd +from flopy.export.shapefile_utils import shape_attr_name from ..datbase import DataInterface, DataListInterface, DataType from ..utils.recarray_utils import create_empty_recarray @@ -124,6 +125,69 @@ def get_empty(self, ncell=0): d = create_empty_recarray(ncell, self.dtype, default_value=-1.0e10) return d + def to_geodataframe( + self, gdf=None, kper=0, sparse=False, truncate_attrs=False, **kwargs + ): + """ + Method to add data to a GeoDataFrame for exporting as a geospatial file + + Parameters + ---------- + gdf : GeoDataFrame + optional GeoDataFrame instance. If GeoDataFrame is None, one will be + constructed from modelgrid information + kper : int + stress period to export + sparse : bool + boolean flag for sparse dataframe construction. Default is False + truncate_attrs : bool + method to truncate attribute names for shapefile restrictions + + Returns + ------- + GeoDataFrame + """ + if self.model is None: + return gdf + else: + modelgrid = self.model.modelgrid + if modelgrid is None: + return gdf + + if gdf is None: + gdf = modelgrid.to_geodataframe() + + data = self.array + + col_names = [] + for name, array4d in data.items(): + if truncate_attrs: + name = shape_attr_name(name, length=4) + else: + name = f"{self.name[0].lower()}_{name}" + array = array4d[kper] + if modelgrid.grid_type == "unstructured": + array = array.ravel() + if truncate_attrs: + aname = f"{name}{kper}" + gdf[aname] = array + col_names.append(aname) + else: + for lay in range(modelgrid.nlay): + arr = array[lay].ravel() + if truncate_attrs: + aname = f"{name}{lay}{kper}" + else: + aname = f"{name}_{lay}_{kper}" + gdf[aname] = arr.ravel() + col_names.append(aname) + + if sparse: + gdf = gdf.dropna(subset=col_names, how="all") + gdf = gdf.dropna(axis="columns", how="all") + + return gdf + def export(self, f, **kwargs): from .. import export