issue import modelgrid from gridgen

[MichelCraninx]

Hello
I have set-up a quadtree model with gridgen (disv-object). Because this takes a long time to calculate I would like to save the modelgrid and import it. For this I have run a simulation so I can import the disv-object and model. When I import the modelgrid the spatial data is missing when I print mg (mg = gwf.modelgrid). However when I plot the modelgrid (mg.plot()), the modelgrid is displayed with the correct coordinates? When I follow the steps: https://flopy.readthedocs.io/en/latest/Notebooks/modelgrid_examples.html, I can change the coordinates to the correct coordinates, however when I plot the new modelgrid the modelgrid is not on the correct position. How can I correct load the modelgrid, resulting from the gridgen calculations?

[martclanor]

Hello. If I understand correctly, this is similar to #2283 which has been fixed since 3.8.0 (https://github.com/modflowpy/flopy/releases/tag/3.8.0). Can you check if this gets fixed upon updating your FloPy to a newer version?

[MichelCraninx]

Hello, I work with version 3.8.2. (recently installed)

[martclanor]

Oh okay, I'm not sure, but can you check if this is also the case before you exported the modelgrid? Just to rule-out that this has something to do with the re-importing of the modelgrid.

I just tested this but on a different grid:

And then setting coord info seems to work fine:

Btw, I am using flopy version 3.8.2 here:

flopy.__version__  # '3.8.2'

[MichelCraninx]

This is how it looks after the import modelgrid from sim:

This is how it looks before importing (original)

The goal is to export the 'base_grid' after gridgen to use for further steps, for example resampling raster data, it works when I use the 'base_grid' directly coming from the gridgen step:

But when I use the modelgrid after importing for further steps: for example resample rasters

then I get this error when I plot the results:


This is not the case when working with the original base_grid ...

It is not clear where the error is, but it seems that importing the modelgrid is causing the error?

[martclanor]

I just realized that your "imported modelgrid" is the one coming from the gwf object which is not aware of the crs because the grid is defined via the ModflowGwfdisv class (which doesn't have the crs parameter).

Does it work if you do something like this upon importing the modelgrid?:

mg_test = gwf_quadtree_test.modelgrid
mg_test.set_coord_info(crs="EPSG:31370")

[MichelCraninx]

Hello, no, I get still the same error, when I only adjust the crs, the plot of the mg_test seems on the wright position, but the print gives = 0, and I cannot do the resampling, with error unsupported operand type, when I both correct coordinates and crs, the plot is on the wrong location and I get the error Raster and modelgrid do not intersect.

[martclanor]

Oh okay. I'm starting to run out of ideas 😅.

I'm wondering why printing the VertexGrid object in both cases show that xoffset (xll) and yoffset (yll) are zero. Did you also pass the correct xorigin and yorigin in your ModflowGwfdisv object?

[jlarsen-usgs]

@martclanor and @MichelCraninx

When gridgen builds a VertexGrid from a StructuredGrid that has offsets applied to it (in FloPy_, gridgen will calculate the vertices for the new VertexGrid from the offset original grid. Therefore, the VertexGrid's vertices are specified in the correct place and no offsets need to be applied to the new grid. Hense, xoff=0, yoff=0, angrot=0 when you load the new modelgrid. If you then apply the offsets to it, you are translating vertices that are in the correct place by whatever offset distance you applied.

In short, xoffset and yoffset should both equal 0 on your new grid (do not apply offsets to it) as the actual vertices that are defined for each cell correspond to the geographic coordinates of the StructuredGrid it was created from.

[MichelCraninx]

Thank you very much martclanor and jlarsen! So I have to find a manner to load correctly the base_grid after gridgen. When I print the properties of the grid after import from the sim object, it is a "vertexgrid" with fine and rough resolutions, as constructed during with gridgen. However, I cannot explain the errors when I use the modelgrid for resampling for example: "unsupported operand type". It is already good to know that is not related with the coordinates.

[jlarsen-usgs]

@MichelCraninx
It's hard to say what's going on without seeing the full error message here. Which specific function call is leading to the ValueError? What's the type() of the rA0220_thickness_clip?

[MichelCraninx]

Hi
I make a modelgrid with gridgen and couple the base_grid to the disv object:

For saving the modelgrid I start a simulation and load the results in a sim:

The modelgrid looks ok:


Then I do a resampling with this modelgrid:


When I plot the resampled plot I get this error:

In case I work directly from the base_grid after gridgen it is ok:

I hope that there is a way I can import the modelgrid after gridgen because we would like to increase the resolution, which will take some time.

Thank you very much!
Michel

[jlarsen-usgs]

I see what's going on. Your base_grid doesn't have the nlay parameter set for it. Here is one way to fix it in your code:

if local_grid_option == 1:
    gridprops = g.get_gridprops_vertexgrid()
    base_grid = flopy.discretization.VertexGrid(nlay=1, **gridprops)

Also rlayer_top is a numpy array that is resampled from your raster. As such, it will not have a head() method. Only flopy.mf6.ModflowGwf and flopy.mf6.ModflowGwfoc objects have the head() method. Example:

gwf = sim.get_model("my model name")
head_obj = gwf.head()
head_array_data = head_obj.get_alldata()  # 4d array of (time, nlay, nrow, ncol)