Feature Request? GridIntersect filter for Active Domain Cells

[javgs-bd]

I'm working with a model construction workflow where it would be great to have a builtin filter for active cells in the Gridintersect.intersect methods. This is, by means of an boolean switch argument be able to exclude from the intersection return inactive cells.
I went through the Flopy docs but seems to me that the filter_query methods available do not enable a direct post-intersection filter of this sort. I wonder if this makes sense or if it is too obvious to do it a posteriori and I'm missing one simple python command.
Any hints would be appreciated.

[dbrakenhoff]

Hi @javgs-bd, thanks for the suggestion. I think that would be a logical feature to add. I'll keep this in mind when I take a look at GridIntersect some time soon (I've been meaning to implement pandas DataFrames for the intersection results).

As for another method, the way I've gone about it thus far, is that I usually end up creating arrays from the intersection results and then using the idomain array (>0) as a mask to update that array.

[javgs-bd]

Hi @dbrakenhoff.

Thanks for the reception.

Yes, for now it takes me only one extraloop across the intersection result.

This is a very rough workaround, as I'm working with a structured grid where idomain is constant through layers. I lost the rest of the recarray fields though. I think there are other numpy methods to drop records based on conditions, but am in a hurry and did not want to do recursive slicing for fields I don't need rigth now..

Just in case somebody else is looking for syntax in the meantime:

active_intersection = []

for row, col in intersection.cellids:
    
    if gwf.dis.idomain.array[0,row,col]==1:
     
        active_intersection.append((row,col))

intersection = np.array(active_intersection)

[rhugman]

Hi @dbrakenhoff.

Thanks for the reception.

Yes, for now it takes me only one extraloop across the intersection result.

This is a very rough workaround, as I'm working with a structured grid where idomain is constant through layers. I lost the rest of the recarray fields though. I think there are other numpy methods to drop records based on conditions, but am in a hurry and did not want to do recursive slicing for fields I don't need rigth now..

Just in case somebody else is looking for syntax in the meantime:

active_intersection = []



for row, col in intersection.cellids:

    

    if gwf.dis.idomain.array[0,row,col]==1:

     

        active_intersection.append((row,col))



intersection = np.array(active_intersection)

I find it easier to pass the intersected rec array to a Pandas dataframe (maybe because im just more used to pandas than recarrays...) and then filter the df with active cellids. Avoids a for loop.

[dbrakenhoff]

Just to add to the above, this is how I sometimes go about it:

gr = GridIntersect(modelgrid)
result = gr.intersect(polygon)

# some array to store intersection cells
mask_intersection = np.zeros(modelgrid.shape, dtype=int)

# structured grid
rows, cols = list(zip(*result.cellids)) # convert to rows and column lists
mask_intersection[:, rows, cols] = 1

# disv grid
mask_intersection[:, result.cellids.astype(int)] = 1

# active intersection
active_intersection = mask_intersection & (idomain > 0)