GridIntersect fails with specific shapely Polygon due to change of nrow, ncol in the discretization object

[javgs-bd]

The code below works if ncol = 50 and nrow = 80, but fails if ncol = 25 and nrow = 40. Am still trying to figure out why, but in any case I though it should be reported.

This is the result of the version printing:

3.8.8 (default, Apr 13 2021, 15:08:03) [MSC v.1916 64 bit (AMD64)]
numpy version: 1.20.1
flopy version: 3.3.5
matplotlib version: 3.3.4
shapely version: 1.7.1

Result with nrow=50, ncol=80

Error with nrow=25, ncol=40

KeyError Traceback (most recent call last)
C:\Users\XXXXX~1\AppData\Local\Temp/ipykernel_20404/762291030.py in
74
75
---> 76 intersection = gridint.intersect(shply_polygon)
77
78

~.conda\envs\flopy-develop\lib\site-packages\flopy\utils\gridintersect.py in intersect(self, shp, **kwargs)
231 and self.mfgrid.grid_type == "structured"
232 ):
--> 233 rec = self._intersect_polygon_structured(shp)
234 else:
235 rec = self._intersect_polygon_shapely(shp, sort_by_cellid)

~.conda\envs\flopy-develop\lib\site-packages\flopy\utils\gridintersect.py in _intersect_polygon_structured(self, shp)
1345 )
1346 else:
-> 1347 v_realworld = intersect.geo_interface[
1348 "coordinates"
1349 ]

KeyError: 'coordinates'

Code

import sys
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as clt
import flopy

from flopy.utils.gridintersect import GridIntersect
from flopy.utils.geometry import shape
import flopy.discretization as fgrid

import shapely
from shapely.geometry import Point, MultiPoint, LineString, Polygon


print(sys.version)
print('numpy version: {}'.format(np.__version__))
print('flopy version: {}'.format(flopy.__version__))
print('matplotlib version: {}'.format(mpl.__version__))
print('shapely version: {}'.format(shapely.__version__))


points  = [
(-18.261405672010028,493.7361282367448),
(20.0,490.0),
(80.0,490.0),
(120.0,450.0),
(150.0,380.0),
(160.0,350.0),
(160.0,310.0),
(160.0,290.0),
(150.0,280.0),
(140.0,250.0),
(130.0,240.0),
(120.0,230.0),
(90.0,220.0),
(90.0,220.0),
(70.0,240.0),
(60.0,290.0),
(30.0,280.0),
(20.0,250.0),
(20.0,250.0),
(20.0,210.0),
(20.0,170.0),
(10.0,160.0),
(10.0,160.0),
(0.0,150.0),
(-22.404438964241834,169.54377311960536),
(-70.04932182490768,284.51294697903825),
(-79.37114673242925,394.30332922318127),
(-83.51418002466107,444.019728729963)
]

shply_polygon = Polygon(shell=points)

shply_polygon


Lx      = 500.              
Ly      = 800.              

nlay    = 1                
ncol    = 50
nrow    = 80           
             
delr    = np.ones(ncol)*Lx / ncol        
delc    = np.ones(nrow)*Ly / nrow    

sgr = fgrid.StructuredGrid(delc, delr, top=None, botm=None, xoff=0, yoff=0, angrot=0)  


gridint  = GridIntersect(sgr, method="structured", rtree=False)


intersection  = gridint.intersect(shply_polygon)



fig, ax = plt.subplots(1, 1, figsize=(8, 8))

sgr.plot(ax=ax)

gridint.plot_polygon(intersection,ax = ax, lw=2) 

[javgs-bd]

I have removed the repeated points in the polygon and still get the error.

[dbrakenhoff]

Thanks for posting this issue. This is a bug in the structured grid intersection routines.

The problem is that in one of the intersection results yields a GeometryCollection containing a Polygon and a Linestring. The shape overlaps the grid cell but then there is also a portion of the polygon that follows the cell boundary for a bit. The only check that is performed is whether the area of the result is larger than 0.0 (which it is). Then it crashes when trying to interpret the resulting GeometryCollection as a Polygon.

I have a fix that I will submit soon, but I may have to wait a bit for #1212 first(?). This fix filters the intersection result to only include Polygons.

For now the solution is to avoid the structured routine, so just use the vertex method:

gridint =  GridIntersect(sgr, method="vertex")
intersection = gridint.intersect(shply_polygon)

[javgs-bd]

Many thanks @dbrakenhoff
Will do that in the meantime.

If I may ask on a parallel topic, that might be related to this behavior. I have noticed that the intersection routines are meant to return any cells touched by the geometry feature. Also, work is being done to extend this behavior in #1198. Which is great, I guess, to achieve things like permeability continuity from line intersections in faults representation.

However, with polygons I believe, an optional, slightly different behavior would be useful. Currently, when dealing with adjacent polygons, intersections will yield shared/repeated cells, since the requisite is only that the intersected area be greater than zero. I believe this can be seen in gridintersect.py lines 1311 to 1315 for the structured case.

If that minimum percentage of area required for a polygon to intersect a cell would be 0.5001, or something representing a bit more than 50%, then the intersection by adjacent polygons would be much more easy to automate.

Is this a reasonable feature to request for the GridIntersection routines related to polygons?

[dbrakenhoff]

I think I understand your question. For adjacent polygons you want the intersection result for each grid cell to have only entry, which could be one of several options:

• the polygon with the largest intersected area (max_area),
• an intersection area larger than 50% of the cell area (min_area),
• or e.g. the one that contains the cell midpoint (midpoint).

Question is whether we want to support that directly, or whether this is a post-processing step. For the intersection() method it kind of makes sense to just return all intersections. And at the moment there isn't a great way to identify polygons within a multi-polygon, which makes it hard to do something similar to the first option directly. However, we could support different methods for storing the intersection result according to the next 2 rules, i.e. min_area = 0.5, or c.contains(midpoint). Something to think about...

As for a postprocessing step, here is a code snippet that does something similar to what you're referring to, using pandas:

# intersect grid with adjacent polygons p1 and p2
# do this separately to avoid getting MultiPolygon intersection results
result1 = gridint.intersect(p1)
result2 = gridint.intersect(p2)

# convert to pandas.DataFrame
df = pd.concat([pd.DataFrame(result1), pd.DataFrame(result2)])

# groupby cellids, and select largest area, or do something else
for cid, group in df.groupby(by="cellids"):
    # do something with largest ...
    largest = group.areas.idxmax()  # (max_area)
    # OR
    # do something with mask that selects entry that covers >50% of cell area 
    cell_area = ... # get cell areas from grid obj
    mask_halfcovered = group.areas > 0.5 * cell_area.loc[cid]  # (min_area=0.5)
    # OR
    # do something with mask that selects polygon that contains midpoint
    xmid = gr.xcellcenters[cid]  # get these from grid obj
    ymid = gr.ycellcenters[cid]  # get these from grid obj
    mask_midpoint = [p.contains(Point(xmid, ymid) for p in group.ixshapes]  # (midpoint)

[javgs-bd]

Thank you very much @dbrakenhoff. You got it perfectly. I was actually thinking on the lines of option 2 and 3. I guess option 3 seems the most direct and easy to understand.
I believe it would be useful not only with adjacent polygons, but also in other instances, like when defining active domains. I've faced this issue with the gridgen intersection methods that behaves similar to GridIntersect. I guess one would not want a cell included in the active domain if the polygon intersection is a tiny fraction of the area.
In anycase, it would be good to have more opinions on this.

[dbrakenhoff]

I'll submit a PR with these additions somewhere this week.