GridIntersect class for unstructured grids

[stub68]

I am trying to intersect an unstructured grid with one or more polygons. I see there is a GridIntersect class in flopy.utils which should do the job nicely (and with a min_area_fraction). However it seems not to have been implemented for unstructured grids (see error below). Is this correct? And is it something that can easily be included? I can achieve the same thing using geopandas sjoin, but need extra steps to select nodes by a minimum area fraction. Intersecting straight from the grid object would be more efficient, but its not a show-stopper. Thanks.

NotImplementedError Traceback (most recent call last)

NotImplementedError Traceback (most recent call last)
Cell In[14], line 2
1 # Create the GridIntersect class for our modelgrid
----> 2 ix = GridIntersect(grid, method = 'vertex')

File ~\anaconda3\lib\site-packages\flopy\utils\gridintersect.py:179, in GridIntersect.init(self, mfgrid, method, rtree, local)
174 self.rtree = rtree
176 # really only necessary for method=='vertex' as structured methods
177 # do not require a full list of shapely geometries, but useful to be
178 # able to obtain the grid shapes nonetheless
--> 179 self._set_method_get_gridshapes()
181 if self.method == "vertex":
182 # build arrays of geoms and cellids
183 self.geoms, self.cellids = self._get_gridshapes()

File ~\anaconda3\lib\site-packages\flopy\utils\gridintersect.py:345, in GridIntersect._set_method_get_gridshapes(self)
343 self._get_gridshapes = self._vtx_grid_to_geoms_cellids
344 elif self.mfgrid.grid_type == "unstructured":
--> 345 raise NotImplementedError()

NotImplementedError:

[christianlangevin]

FloPy supports three different types of model grid objects to varying degrees. These grids are StructuredGrid, VertexGrid, and UnstructuredGrid. These model grids are typically created under-the-hood by flopy when you create a model. StructuredGrid corresponds to a DIS-type MODFLOW grid and is well supported throughout flopy. VertexGrid corresponds to a DISV-type MODFLOW grid in which the mesh is the same for all layers. A VertexGrid should be supported by the grid intersect classes. The DISU-type UnstructuredGrid object is probably the least well supported, though support has been improving. I suspect you are using the DISU UnstructuredGrid. If so, you might be able to create a separate VertexGrid object for your model (by instantiating a VertexGrid object manually) and use it for the intersects. We plan to continue to add better support for DISU unstructured grids, but the flexibility to change the mesh by layer, or not use layers at all, poses some challenges for full support.

[stub68]

Thanks Chris. Yes it's a DISU grid for a USG model. Im concerned that node numbering order may differ between DISU and DISV. But in any case, I should be able to do the intersection using a geodataframe. Thanks again.

[ccarson-tt]

I used a geopandas workaround for this as well - process was like this:

1. Read in shapefile with geopandas
2. Create df containing coordinates for all nodes in modelgrid
3. Use Polygon.contains(point) to determine if the node is within in polygon bounds - also used another small function to assign cells to zones in my model grid that weren't within any of the polygons (find_closest_polygon).

Code:

import geopandas as gpd
from shapely.geometry import Point
import pandas as pd

def find_closest_polygon(point, polygons):
distances = [polygon.distance(point) for polygon in polygons]
closest_index = np.argmin(distances)
return closest_index

zone_list = []

for index, row in xy_df.iterrows():
point = Point(row['x'], row['y'])
found_polygon = False
for ipolygon in range(len(zones.geometry)):
if zones.geometry[ipolygon].contains(point):
zone_list.append(ipolygon+1)
found_polygon = True
break
if not found_polygon:
closest_index = find_closest_polygon(point, zones.geometry)
zone_list.append(closest_index + 1)

Cheers,

Colter

[jlarsen-usgs]

Just as a follow up. The grid objects now have a geo_dataframe attribute that will return a geopandas dataframe of the geometries. The syntax is:

modelgrid = mymodel.modelgrid
gdf = modelgrid.geo_dataframe

[dbrakenhoff]

Adding DISU support to GridIntersect requires a list of geometries and node numbers from the DISU grid, which isn't that difficult to implement. The challenging part is how to deal with the 3rd dimension. Some questions/thoughts to answer to figure out how we could support this:

• Do we add all DISU cells at all elevations for intersections? Then the user has to post-process the intersection results to obtain the cells at the correct elevations. This seems fairly simple to implement.
• Or can we add z data (e.g. elevation of cell center) to the geometries (that are generally ignored in geometric operations as far as I know)? That z-data is persisted to the intersection results, so maybe that information can maybe be used to filter the result.
• If a DISU maintains some kind of layering, somehow determine a mapping between layer number and node numbers? That mapping can be used to select the correct cells for intersection.
• Add some methods to help users filter results based on elevation constraints?

Maybe worth starting an issue to discuss this further?

[stub68]

Thanks for the replies and suggestions. @ccarson-tt: Yes I have used a similar workaround using geopandas, with and additional few lines to select based on the minimum overlap area for each cell/node. @jlarsen-usgs - Interestingly this approach is faster than using the .geo_dataframe attribute in flopy, at least when I tried it.

@dbrakenhoff: Regarding the possible implementation of GridIntersect for DISU, it would be useful to add the z elevation to the resulting geodataframe, since, as you say, further post processing would be needed for the layered cells/nodes. I'm not sure that DISU includes the concept of layers. If it does, the a layer attribute would also be very useful.