Enable a mechanism to capture coordinate reference system metadata

[mwtoews]

This is a feature request to enable a MF6 model to capture information for a coordinate references system (CRS; also called spatial reference system or SRS), so this information can be used by GIS software and tools. This information would not be used directly by MODFLOW 6, but could potentially be passed to the binary grid file (.grb).

Background

There are several mechanisms that software can capture CRS info. The simplest method is to use look-up authority names and codes. The most commonly used authority is EPSG, but several authorities also exist. For example, using pyproj:

$ python -c "import pyproj; print(pyproj.database.get_authorities())"
['EPSG', 'ESRI', 'IAU_2015', 'IGNF', 'NKG', 'OGC', 'PROJ']

and the second piece is the code, which is always an integer greater than one. For example, here are some projected CRSs on Earth:

• EPSG:26916 - NAD83 / UTM zone 16N
• ESRI:102007 - NAD 1983 Albers Hawaii
• IAU:2015:39916 - Earth (2015) / Ographic/ Equirectangular, clon = 180

Software like PyPROJ can identify these CRS like this:

from pyproj import CRS
crs1 = CRS.from_authority("EPSG", 26916)
crs2 = CRS.from_authority("ESRI", 102007)
crs3 = CRS.from_authority("IAU_2015", 39916)

and Esri's proprietary ArcPy has a SpatialReference interface:

import arcpy
sr1 = arcpy.SpatialReference(26916)
sr2 = arcpy.SpatialReference(102007)

(according to their docs, the WKID code will either refer to EPSG or ESRI, but there is no need to specify the authority since the codes don't overlap.)

New authority names and codes are added/updated every year, so capturing this data is only as good as the client software is able to look up the information (e.g. older software won't understand future codes).

An alternative mechanism to capture CRS is to explicitly capture the coordinate reference types and parameters, which is complicated (e.g. WKT1, WKT2:2019, PROJJSON, CF Metadata Conventions, etc.). This complex approach should be avoided due to burden of information that would need to be handled, and the complexity of different markup schemes.

This approach would limit users to only use registered CRS codes. Custom or local CRSs would not be permitted.

Discussion of idea

A suggestion for discritization input option blocks is to be modified as follows:

BEGIN OPTIONS
  [LENGTH_UNITS <length_units>]
  [NOGRB]
  [CRS_DEF <crs_def>]
  [XORIGIN <xorigin>]
  [YORIGIN <yorigin>]
  [ANGROT <ang>]
END OPTIONS

This would describe the short CRS definition for XORIGIN/YORIGIN, but have no other relation with LENGTH_UNITS or ANGROT, which are local coordinates used by MODFLOW 6.

The new option CRS_DEF would be a string with some valid entry examples:

EPSG:26916
ESRI:102007
IAU:2015:39916
FUTURE:123

invalid entries would be missing an authority name (e.g. 26916), missing ":" separator (e.g. EPSG_26916), or have non-integer authority code (e.g. EPSG:26916A).

Using some Python pseudocode, this data could be stored as two entities:

crs_def = "IAU:2015:39916"
crs_auth_name, crs_auth_code = crs_def.rsplit(":", 1)
crs_auth_code = int(crs_auth_code)
print(f"{crs_auth_name=}, {crs_auth_code=}")
# crs_auth_name='IAU:2015', crs_auth_code=39916

within Fortran, crs_auth_name could be a character string with max length 10 and crs_auth_code would be an integer. Or crs_def could be stored as a single character string with max length 20.

Unsure if/how this would make it into the binary grid specification. The current is "VERSION 1", so if a next version is planned, this information could be stored in the .grb file. I've also noticed that the format specification is missing LENGTH_UNITS too.

[langevin-usgs]

@mjreno, drawing your attention to this issue as @mwtoews has thoughts on how to specify the model CRS. Not sure this is consistent with your NetCDF work.