How to get l,r,c of a cell that contains coordinates x, y and z using flopy.mf6

[chenmingzhang]

This function is very handy during mesh dependency test, as we do not have to change the lrc of the wells during remeshing.

I used to use the follow function to achieve this, if mf2005 is used :

dis = flopy.modflow.ModflowDis( ....)
def get_lrc_from_coordinates(x,y,z,dis=dis) :
    """
    function to get the lrc from the coordinates. 
    """
    [row,column]= dis.get_rc_from_node_coordinates(x,y)
    layer = dis.get_layer(row,column,z)
    return (layer, row , column)

How to achieve this if using mf6, given dis.get_rc_from_node_coordinates and dis.get_layer are not available?

Following the examples, I have done:

from flopy.utils.gridintersect import GridIntersect
from shapely.geometry import Point
ix = GridIntersect(gwf.modelgrid, method='vertex')
rc_point=ix.intersects(Point(150,150))
l_point = ??????? (rc_point, z= -10)   

where ????? is needed. or my way is completely wrong.

I have checked the methods in flopy.mf6.ModflowGwf.modelgrid and it seems nothing is for this.

Thanks, Chenming

[spaulins-usgs]

@chenmingzhang, did you try the modelgrid's intersect method? You can do try something like:

layer, row, column = gwf.modelgrid.intersect(x, y, z, local=True)

If you want to define a coordinate system you will first want to set the coordinates with the modelgrid's "set_coord_info" method and then call the intersect method, setting local to False (default). Note that intersect returns different results for vertex and unstructured grids.

[chenmingzhang]

@spaulins-usgs thanks for the reply. Following your advice, I tried a few as below:

1:

layer, row, column = gwf.modelgrid.intersect(x, y, z, local=True)
--------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [70], in <cell line: 1>()
----> 1 layer, row, column = gwf.modelgrid.intersect(x, y, z, local=True)

TypeError: intersect() got multiple values for argument 'local'

2:

layer, row, column = gwf.modelgrid.intersect(x, y, z)
print (layer,row,column)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [73], in <cell line: 1>()
----> 1 layer, row, column = gwf.modelgrid.intersect(x, y, z)
      2 print (layer,row,column)

ValueError: not enough values to unpack (expected 3, got 2)

3:

lrc = gwf.modelgrid.intersect(x, y, z)
print (lrc)
(1, 5)

Which workds while not provide z.

I had a looked into the code and saw:

File ~/python-3.9/lib/python3.9/site-packages/flopy/discretization/grid.py:673, in Grid.intersect(self, x, y, local, forgive)
    671 def intersect(self, x, y, local=False, forgive=False):
    672     if not local:
--> 673         return self.get_local_coords(x, y)
    674     else:
    675         return x, y

suggesting that z is not a input argument for fuction intersect.

I am aware the discussion , which seems not in flopy 3.3.6 that i am using.

Cheers, Chenming

[jlarsen-usgs]

@chenmingzhang

Update your flopy installation to the most recent development version and the code snippet that spaulins-usgs shared will work

pip install -U https://github.com/modflowpy/flopy/zipball/develop