Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Vertex coordinates are not relative to the specified origin #2388

Open
cneyens opened this issue Dec 2, 2024 · 0 comments
Open

Vertex coordinates are not relative to the specified origin #2388

cneyens opened this issue Dec 2, 2024 · 0 comments
Labels

Comments

@cneyens
Copy link

cneyens commented Dec 2, 2024

When using a Vertex grid object, the values of the vertex coordinates (verts instance) change based on the specified origin. They are given in georeferenced coordinates, i.e. the sum of the origin and the vertex coordinates in the model projection + optional rotation. I find this confusing, as this differs from what is written to the DISV file. Is this intended behaviour? I would expect that the vertex coordinates are always relative to the model origin.

For example, when reading in a model with a specified origin, the vertex coordinates in the modelgrid object are already offset by the origin and do not resemble the values in the DISV file on disk. When re-using those values, e.g. to create a DISV grid for a different model, the origin should not be taken into account anymore as that would double the offset (see example below). But then of course the origin information is no longer contained in the new DISV file.

Possibly related to #2283; see also #2364.

Example code
import os
import numpy as np
import flopy
from flopy.utils.triangle import Triangle

# create triangular grid
triangle_ws = 'triangle'
if not os.path.exists(triangle_ws):
  os.mkdir(triangle_ws)

active_area = [(0, 0), (0, 1000), (1000, 1000), (1000, 0)]
tri = Triangle(model_ws=triangle_ws, angle=30)
tri.add_polygon(active_area)
tri.add_region((1, 1), maximum_area=50**2)

tri.build()

# build vertex grid object
vgrid = flopy.discretization.VertexGrid(
  vertices=tri.get_vertices(),
  cell2d=tri.get_cell2d(),
  xoff=5555,
  yoff=7777,
  crs=31370,
  angrot=0,
)

# create MODFLOW 6 model
sim = flopy.mf6.MFSimulation(sim_name='prj-test', sim_ws='model')
tdis = flopy.mf6.ModflowTdis(sim)
ims = flopy.mf6.ModflowIms(sim)

gwf = flopy.mf6.ModflowGwf(sim, modelname='gwf')
disv = flopy.mf6.ModflowGwfdisv(
  gwf,
  xorigin=vgrid.xoffset, # specified origin
  yorigin=vgrid.yoffset,
  angrot=vgrid.angrot,
  nlay=1,
  top=0.0,
  botm=-10.0,
  ncpl=vgrid.ncpl,
  nvert=vgrid.nvert,
  cell2d=vgrid.cell2d,
  vertices=tri.get_vertices() 
)

print(disv.vertices.get_data()[:5]) # no offset
print(gwf.modelgrid.verts[:5]) # includes offset 
print(gwf.modelgrid) # origin is set in modelgrid, so verts is no longer relative to this origin
disv.write() # no offset in written vertices coordinates

# creating a new DISV object from this modelgrid is confusing
# let's assume I loaded the gwf-model from disk, so I don't have access to the the Triangle output,
# only the gwf.modelgrid object
grid = gwf.modelgrid

prt = flopy.mf6.ModflowPrt(sim, modelname='prt')
disvprt = flopy.mf6.ModflowPrtdisv(prt,
                                 xorigin=grid.xoffset, # specified origin
                                 yorigin=grid.yoffset,
                                 angrot=grid.angrot,
                                 nlay=grid.nlay,
                                 top=grid.top,
                                 botm=grid.botm,
                                 ncpl=grid.ncpl,
                                 nvert=grid.nvert,
                                 cell2d=grid.cell2d,
                                 vertices=np.column_stack((np.arange(grid.nvert), grid.verts)).tolist()
                                 )
print(prt.modelgrid) # offset
# Doubled offset in verts. I can remove the origin in the creation of DISV above, but then it would not write this info to disk.
# With a DIS grid, I do need to specify the origin coordinates.
print(prt.modelgrid.verts[:5]) 

Expected behavior
I would expect the modelgrid verts instance to remain static when called by the end-user and always resemble the coordinate values relative to the specified offset, as is written to the DISV file. Under the hood, this may change for e.g. plotting. Alternatively, the offset in the modelgrid object should change depending on whether or not the verts instance resembles relative or absolute coordinates, or an extra instance with the static relative vertex coordinates could be added to the modelgrid object (perhaps that's the easiest option?).

Flopy v3.8.2

@cneyens cneyens added the bug label Dec 2, 2024
@cneyens cneyens changed the title bug: Vertex coordinates are not relative to the specified origin Vertex coordinates are not relative to the specified origin Dec 2, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

1 participant