How does Flopy locate models now in reference coordinates when the computational units are in feet?

consider this model, which is in feet, but referenced to a coordinate reference system that is in meters:

import os
import numpy as np
import flopy
fm = flopy.modflow
from flopy.discretization import StructuredGrid
from flopy.utils.reference import SpatialReference

dis_kwargs = {#'nlay': 1,
              'nrow': 10,
              'ncol': 10,
              'delr': 1000, # in feet
              'delc': 1000 # in feet
}

# grid location info
xoff = 1000 # lower left x-coordinate in UTM meters
yoff = 1000 # lower left y-coordinate
rotation: 0.
epsg = 26916

tmpdir = 'tmp'
if not os.path.isdir(tmpdir):
    os.makedirs(tmpdir)
    
mf2005 = fm.Modflow('mf2005.nam', model_ws=tmpdir)
dis = fm.ModflowDis(mf2005, **dis_kwargs, lenuni=1)

Previously, with SpatialReference, one could just specify the lower left or upper left corner, and an epsg code for a CRS, and it would convert the grid to the correct location in the CRS:

mf2005.sr.xll = xoff
mf2005.sr.yll = yoff
mf2005.sr.epsg = epsg
mf2005.sr.bounds
(1000.0, 1000.0, 4048.0, 4048.0)

Unless I'm missing something, with modelgrid, there doesn't seem to be a way to do this:

mf2005.modelgrid.set_coord_info(xoff=xoff, yoff=yoff, epsg=epsg)
mf2005.modelgrid.extent
(1000.0, 11000.0, 1000.0, 11000.0)

mg = StructuredGrid(delc=np.ones(10) * 304.8, 
                                    delr=np.ones(10) * 304.8, 
                                    xoff=xoff, yoff=yoff, angrot=0, 
                                    epsg=26916, lenuni=2)
mf2005.modelgrid = mg
mf2005.modelgrid.extent
(1000.0, 11000.0, 1000.0, 11000.0)

(the lenuni argument doesn't seem to have any effect on the last result)

There also doesn't seem to be a way to supply the upper left corner instead of the lower left corner. I get that this makes things a little more complicated, but oftentimes (as in the USGS archiving guidence), only the upper left is given.

SpatialReference supported this as well

sr = SpatialReference(delr=mf2005.dis.delr, 
                      delc=mf2005.dis.delc,
                      xul=1000, yul=4048, epsg=26916, 
                      lenuni=1 # model lenuni
)
sr.yll
1000.0

there is an an _yul_to_yll() method, but it only works if the CRS and model units are the same:

mf2005.modelgrid._yul_to_yll(4048)
-5952.0

[christianlangevin]

You're right, @aleaf-usgs, this is not possible now with modelgrid, and it's by design. But I also ran into this same issue this past week. The thinking is that if you are building a model in feet that you should be using a projection that is also in feet. But the user that I was working didn't feel like that was a fair requirement, and I can see that point. What I did in that case was create another another model grid object (a StructuredGrid object) but I specified my DELR and DELC values to be in the real-world map units. I then used that modelgrid object for the plotting. I personally don't like the LENUNI value being used to make that conversion automatically because it seems too hidden. After working with this stuff a bit the past week, I could how it might be beneficial to allow the user to swap out the modelgrid object with a custom one that would be used for all plotting and exporting, but we wouldn't want that modelgrid object to be for anything where a length unit might be returned for model processing (intersect, for example).

@langevin-usgs I know you guys have thought about this, and I get the desire to avoid complexity, but as it is now, none of the export stuff (shapefiles or rasters) will work for anyone with a model in feet. Using a second model grid object for the built in export methods doesn't work because they use the attached model grid object.

All the more reason to use meters for new models, but the vast majority of existing models I've had to work with (and I suspect this is true for many in the U.S.) are in feet. Conversely, commonly used projected CRS (UTM, Albers, state coordinate systems, etc.) are usually in meters, even in the U.S. Working with projected CRS that are in feet is generally a huge headache and makes it difficult to share data with people.

In the interest of simplifying things, I think it would make more sense to have a representation of the grid that is always in meters (the global standard for projected CRS), have a default length conversion of 1, and then make it clear in the documentation that if the computational units aren't in meters, it is the user's responsibility to supply the model length units (using lenuni or more human-readable text). Users will pick up on this right away when their grid doesn't show up in the right place, and look for a way to fix it.

[christianlangevin]

@aleaf-usgs, I understand what you are saying. Working with feet-unit models in the U.S. poses some additional challenges. If we were allow for another CRS meter representation, then it would make sense to use that for any map that is constructed or for any export that contains map coordinates. But what about some of the other operations that are less clear? Do we use that meter grid CRS representation for an intersect so that if I pass in a point, it must be in meters? What about for a cross section? Should the z-axis be in feet and the x-distance in meters? Are there other implications that we should discuss? For now, I think all the SpatialReference capabilities are still supported. So the capabilities are still there I think. Perhaps we should shoot for a resolution to this for the following release, since it might be too soon to get this resolved this week for our planned release. We should discuss next week in Golden. Interested flopy users should feel free to comment here if you have an opinion.

@langevin-usgs to answer your questions:

intersections
Yes, I think the default for intersects should be meters.

• Most likely any data that has to go into the model (including field data collected with GPS, data received from other agencies, or online datasets, etc.) is going to either be in a geographic coordinate system and have to be converted to projected coordinates via pyproj or another utility, or be in UTM or a local coordinate system (e.g. WI Transverse Mercator) that is in meters.
• In my experience it is way easier to convert data from GCS to a projected coordinate system if the latter is commonly used (for example, has an epsg code that is recoginized by the PROJ library); the vast majority of these are in meters. This is especially true for users that aren't as experienced with GIS.

cross sections
Ideally, I think this depends on the use case. A site-scale cross section would be more appropriate in the model units. For troubleshooting issues with the grid or locating features in the model, one might want to display row or column for a structured grid. Miles or kilometers might be more appropriate for a report figure documenting a regional model. I think the user should be able to decide what horizontal length units they want to see. It adds complexity, but shouldn't be that hard to implement with a robust unit conversion module and a few limitations on input.

I don't think that modelgrid specifically has to do this if it conflicts with other goals for that class. But I do think that Flopy would be more useful to those of us stuck using english length units if it just "worked" with a minimum of user intervention/headache.

[jlarsen-usgs]

I am working with a user that is having the same issue currently (model in feet, and projection is in meters) who is trying to export shapefiles of model packages. Currently there is no way to override the attached modelgrid class in the model object.

My thought, as a temporary fix, is to either add a method to the shapefile_utils module that would allow the user to export package shapefiles with a modelgrid that they have defined. The call signature would look something like:

def package_to_shapefile(filename, pak, modelgrid):
    """User method to export package shapefiles with a custom modelgrid"""

Does anyone have other ideas/interventions that would make this conversion easy for the user?

[spaulins-usgs]

The shapefile export code supports kwargs. Perhaps the easiest solution is to have model_attributes_to_shapefile (in shapefile_utils.py) use a modelgrid passed in kwargs if it is provided. It already gets epsg and prj strings from kwargs so this would be consistent with that behavior.