Revisions to SpatialReference; generic ModelGrid class

[aleaf]

TLDR: if care about SpatialReference or the way Flopy handles importing, exporting and plotting of geographic information, please comment and provide input! Specifically, input on the purpose and design paradigms of SpatialReference, or any alternative implementations that would better address the issues below if you don't agree with what's proposed.

@langevin-usgs and @jdhughes-usgs have tasked @spaulins-usgs, @jlarsen-usgs and me with refactoring flopy's plotting, export and spatial reference capabilities so that they behave consistently between the MODFLOW-6 and pre-MODFLOW-6 classes, and are also unified/extended to handle unstructured grids. So far, development towards those goals has centered around a generic ModelGrid class, with subclasses for unstructured and structured grids. The ModelGrid family of classes are supposed to provide information about the model grid that is independent of the MODFLOW version or grid input specification (e.g., the cells of any grid can be described by exterior vertices).

One key issue related to both export and plotting is how transformations are handled between model computational coordinates and geographic or "real-world" coordinates. Currently, this functionality is performed for structured grids by the SpatialReference class. The notebook example referenced below provides a brief rundown of the core-functionality and motivation behind SpatialReference. While @jtwhite79 and I (and perhaps others) have found SR to be very useful for pre- and post-processing models, there have also been a number of criticisms including:

• too many attributes
• attribute names are confusing; not clear which ones contain model coordinates vs. geographic coordinates
• no one knows how to use it
• too specific to structured grids
• methods such as export_array() should be in the export module instead

One proposed approach (that is somewhat implemented on the main branch tip for this effort) is to move much of SR's functionality into the ModelGrid class, and have SR only perform transformations. Some problems with this approach are that it doesn't clear up the attribute name confusion, requires more user input (instantiation of both SR class and ModelGrid classes for stand-alone functionality), and it is not backwards-compatible.

Another approach is demonstrated in this Notebook, which also attempts to articulate the main ideas behind SpatialReference. In this revision, the new ModelGrid class stores coordinate information in the model computational coordinates, and SpatialReference transforms and stores the same key attributes in geographic coordinates. Separating the model and geographic coordinate information into two classes with the same attribute interface will hopefully alleviate some of the attribute name confusion.

Note that in response to criticism of attribute names, I am proposing some name changes (which may or may not be helpful):

old name	new name	what it is
xcentergrid	xcenters	array of x-coordinates for model cell centers
ycentergrid	ycenters	array of y-coordinates for model cell centers
xgrid	xedges	array of x-coordinates for model cell corners
ygrid	yedges	array of y-coordinates for model cell corners
get_gridlines()	gridlines	list of tuple sequences that describe cell edges (for plotting)

If we were to change any names, I think at least a 1 release deprecation warning (similar to other python packages) would be the best way to go.

As shown in the above notebook, I am also proposing that we retain export functionality in SR (such as export_array), and even expand SR to including methods for sampling geographic information to the model grid.

Please comment if you have an opinion.

@mwtoews, @mnfienen.

[jtwhite79]

Thanks @aleaf . I can live with the slight attribute name revisions but I suggest we leave the deprecation warnings (and associated backward compat) at least xcentergrid and ycentergrid for a while to allow time for everyone upgrade.

RE that notebook you linked, I think it looks great but I would caution against adding any additional GIS deps to flopy - rather, just expose a clean API for users to access flopy data in various spatial containers (shapefile, raster, netcdf, etc). This is the plan, right?

I still like the option to store the spatial reference info in the existing model input files just so it doesn't require additional files. I think we could easily read this info without requiring a full model load.

Just my thoughts...

[mwtoews]

Thanks @aleaf , finally had a chance to go thru SpatialReference_demo.ipynb and have a look at flopy.grid

• sample(): super handy (the inspiration in Revise SpatialReference properties, add geotransform #367 was to provide geotransform affine transformation coefficients, which is handy for this purpose), but a good implementation may out-weight the lightweight aims described above. If it's a regular grid, then the gdal/rasterio bindings will work for most methods (nearest, bilinear, cubic, average, median, mode, etc.). But for non-uniform grids, gdal/rasterio cannot be used, and code may need to be written to re-implement some of these methods.
• Looking at intersect(), the word "intersect" has a particular meaning for both set theory and GIS, which are different than how it's used in the demo. A better name for this functionality is join() (or sjoin()). But like sample(), this method could get really complex, e.g. searching for the closest to the cell centroid or cell exterior. Or with multiple joins per feature, do we aggregate the average, median, mode, etc. of the joined attributes. There are lot's of different ways a user may want to combine spatial data.
• I'd prefer moving export_array to an export module, with a deprecation warning to the current functionality. And also remove the spline interpolation for Arc Ascii exports with rotated grids.
• If we're renaming things, perhaps move from proj4_str to proj_str. As of the past year, PROJ.4 refers to the project, while PROJ is the library i.e. versions beyond 4. (It was at version 4 for 25 years that everyone seems to have assumed that the version number was part of the name...)
• As for attribute/parameter names that work on local Cartesian coordinate reference systems and projected coordinate references systems,^[1] the current situation isn't too bad. For example, usually "x" means along Eastings for projected, or along rows for local Cartesian, and that seems fine. transform() interchanges the meaning of x based on inverse=True or False, and that's also fine.

^[1] I'm avoiding using "geographic coordinates", as this term often implies degrees of latitude/longitude; length units in degrees are unsuitable for modflow models.

[spaulins-usgs]

@aleaf I apologize for not responding sooner. I have been meaning to reply to your thoughtful description of possible revisions to SpatialReference, but have been out for much of August. My primary concern with SpatialReference is trying to figure out how it can work with MF6's different model grid types. My views on SpatialReference therefore are likely to differ from those who have extensively used flopy's SpatialReference and do not want to see its backward compatibility broken. That said, I think there is good code in SpatialReference and most of what I propose has to do with moving and refactoring that code to work seamlessly with MF6 models (which can have multiple grids and multiple grid types). I understand concerns about backward compatibility and understand that a reasonable case must be made to justify code refactoring.

In early August @langevin-usgs, @jdhughes-usgs, @aleaf, @jlarsen-usgs, and I discussed the SpatialReference issue. While there was not a consensus, there was some interest in having a more bare bones SpatialReference class that primarily does what its name implies. Here is my argument for that and my proposed implementation (I have already implemented much of it).

The problems I see with the current SpatialReference...

1. Currently the primary purpose of SpatialReference is counterintuitive to what it's name implies. This increases the learning curve for new users. Much of SpatialReference currently involves working with a model grid. In @aleaf's description of the purpose of SpatialReference, his first bullet point is "lightweight (doesn't take minutes to load) representation of model grid in geographic coordinates for relating the model to external data". The name SpatialReference does not imply anything to do with a model grid, it implies storing geographic reference information (what you find in the proj4 string). Additionally, our current implementation of SpatialReference is inconsistent with how nearly every other python library implements SpatialReference. Arcpy and other popular libraries have SpatialReference classes that are primarily used to store geographic reference information. The closest analogy I can think of to MODFLOW's model grid is ArcGIS's shapefile. Arcpy's SpatialReference doesn't contain any properties/methods related to shapefiles. Our SpatialReference should not contain any properties/methods related to model grids.

2. Currently parts of the SpatialReference code are implemented to specifically support the structured model grid. For example, optional parameters in SpatialReference's init include delr and delc, which are specific to a structured grid. These variables are pervasive throughout the SpatialReference code, I found 32 instances of delc in SpatialReference. If we continue to support this, I assume we will need to add more optional parameters to support the other grid types as well. This will get very messy very fast. Note that there are other places where SpatialReference assumes a structured grid, like the from_gridspec method. MF6 has three grid types, and there may be more in the future. As support for different grid types is added to SpatialReference, it will become messier and more counterintuitive.

Here is how to propose to solve this...

I implemented a SpatialReference class that primarily contains the proj4, epsg, and crs attributes, and methods specific to these attributes. Most other attributes and methods were moved elsewhere, with much of it moved to the ModelGrid classes described below.

A ModelGrid base class was implemented with initially two subclasses, StructuredModelGrid and VertexModelGrid. Anything specific to a model's grid but not specific to the grid type was moved from the old SpatialReference to the ModelGrid base class, including origin_loc, origin_x, origin_y, lenuni, rotation, and transform (user can get xll, yll, xul, yul from origin info).

The StructuredModelGrid subclass contains attributes/methods from the old SpatialReference class that are specific to the structured model grid including, delc, delr, from_gridspec, getting cell centers, cell vertexes, and cell edges. It also contains other attributes/methods specific to a structured grid including top, botm, and idomain. Note that there is no longer a need to store delc and delr in both SpatialReference and StructuredModelGrid. Instead the information is only in the class with a name that intuitively implies that it belongs there.

The VertexModelGrid subclass is similar to the StructuredModelGrid subclass, but contains attributes specific to that model grid and methods such as getting cell centers, cell vertexes, and cell edges.

The ModelGrid classes can be constructed with or without geographic coordinate information (spatial reference, rotation, origin). If geographic coordinate information is included in the model grid class (passed in in the constructor), the ModelGrid class uses that coordinate system, otherwise it uses model coordinates. I also added an attribute flag that lets you turn off using coordinate reference information even when it is included. If you want information related to the model grid (cell centers, edges, grid lines, etc) you always get it from the model grid class. I think this is more intuitive than getting model grid information from both the model grid class and the spatial reference class.

Note that in many ways the interface of this implementation is very similar to what @aleaf is proposing. I am mostly just proposing that SpatialReference code specific to the model grid is moved from SpatialReference into the appropriate model grid classes. I put together a notebook that mimics the first part of Andy's notebook, which shows you just how similar everything is.

Advantages to this implementation...

• Classes do what their names imply
  -- SpatialReference stores reference information about the spatial coordinate system
  -- ModelGrid base class stores general information about the model grid, how it is rotated, the origin point that it is rotated around, and the units of the grid
  -- StructuredModelGrid and VertexModelGrid contain information specific to those types of model grids
• Functionality is modularized into classes that have very specific well-defined purposes consistent with their names.
• SpatialReference class behaves like the SpatialReference classes in most other popular python libraries. Users familiar with these other popular libraries will have a minimal learning curve.
• Adding new types of model grids will not lead to additional code clutter, just create a new subclass of ModelGrid and put the new code in that subclass
• You can still create a light weight StructuredModelGrid (leave out the optional top, botm, and idomain information) that has the same functionality as Andy's SpatialReference

[jtwhite79]

I tend to agree with @spaulins-usgs that the current SpatialReference is really not extensible in its current form. Its full of misnomers and is not setup to handle unstructured grids. So, it needs to change.

I'm all for the purposed changes, but I (selfishly) need a little time to make the transition. So, is it possible to have wrapper methods in the new, stripped down SpatialReference that trigger a warning but then map the new location of the methods in ModeGrid instance (at least for a release or two)? I'm primarily interested in xcentergrid and ycentergrid in the current SR...

[spaulins-usgs]

@jtwhite79 we can change SpatialReference in two steps. The first step we leave all of the functionality in SpatialReference but add deprecation warnings. Step two we remove the deprecated functionality.

[spaulins-usgs]

Josh and I have implemented some proposed changes to exporting, plotting, spatialreference, and model grids that many of us have already discussed. We are planning to pull into develop in the next two weeks. If you have any comments please provide them in the next week or so.

The implementation is available here:

https://github.com/modflowpy/flopy/tree/spaulins-develop-grid

Exporting and plotting are now implemented using the grid class instead of spatial reference and dis. Both exporting and plotting work the same way on mf2005 and mf6. Backward compatibility has been maintained for both exporting and plotting. Exporting and plotting now also support both structured and vertex grids, with plotting having limited support for unstructured grids. The exporting interface is mostly unchanged (grid object swapped for a spatial reference object), while a new simplified plotting interface has been implemented (the old plotting interface is also still available).

Please take a look at the following notebooks.

These notebooks contain examples of plotting using mf2005 and mf6 using structured, vertex, and unstructured grids. These notebooks also demonstrate the new and old plotting interfaces and show backward compatibility: flopy3_PlotMapView_demo, flopy3_PlotCrossSection_demo

These old notebooks have been updated to show that the plotting code still functions similarly to before: flopy3_MapExample, flopy3_UnstructuredGridPlotting, flopy3_PlotArrayExample.

The exporting interface is mostly unchanged. Examples of exporting with mf2005 can be found in flopy3_export and flopy3_shapefile_export notebooks and exporting with mf6 (same interface) can be found in the flopy3_mf6_tutorial notebook.

The notebooks are all in the github branch here:

https://github.com/modflowpy/flopy/tree/spaulins-develop-grid/examples/Notebooks

There are a few implementation details worth mentioning. For plotting ModelMap has been deprecated and replaced with MapView. When you call the deprecated ModelMap or ModelCrossSection class it gives you a deprecated MapView/CrossSection object back with the old deprecated interface. When you call PlotMapView, you get the appropriate MapView object (Structured, Vertex, or Unstructured) with the new interface. By implementing plotting this way, backward compatibility is maintained alongside a new simplified interface without the need for a plot6 folder. The biggest difference, for the user, between the old interface and the new one is a reduction in the number of optional parameters that can be supplied to instantiate the plotting classes (PlotMapView and PlotCrossSection). The optional parameters xul, yul, xll, yll, rotation, and length_multiplier have been removed from the new interface. However, these parameters can still be supplied to the depreciated ModelMap and ModelCrossSection classes as before.

The exporting utility methods (in shapefile_utils) have the same interface except sr (spatialreference) is replaced with mg (grid). If the code detects a spatialreference object passed in the place of the grid, it gives a deprecation warning, but still works.

Model, package, and data interface base classes were added that both the mf2005 and mf6 model/package/list/array classes derive from. ModelInterface explicitly defines a common interface supported by all types of models (mf2005, mf6, ...). PackageInterface explicitly defines a common interface for all package classes of all types of models. Similarly DataInterface explicitly defines a common interface for all types of data across all types of models. Having these common interfaces along with the new grid classes allowed us to implement plotting and exporting code that is shared between mf2005 and mf6.

[jtwhite79]

@spaulins-usgs - this all looks good to me - until I start actually digging in with use cases, I can't really say much. I do think we need to have an explicit transition plan about when the deprecation warnings become errors (again being selfish because this will break things in pyemu and I suspect will break things for other codes that rely on flopy) - it would be nice to coordinate that change from warnings to errors.

One small thing: why is BaseModel's exe_name attribute now (semi-)protected with the leading underscore? It seems to me we want to leave this available to be changed at will.

[spaulins-usgs]

@jtwhite79 Yes we should definitely coordinate the change from warnings to errors. Concerning exe_name, I see your point. exe_name has now been changed back to its original unprotected state.