Coordinate Reference Systems

Every Survey must have a coordinate reference system (CRS) defined and all datasets within the Survey adhere to the same CRS.

This example explores the CRS variable and shows how it is linked to data variables.

Dataset Reference: Minsley, B.J, Bloss, B.R., Hart, D.J., Fitzpatrick, W., Muldoon, M.A., Stewart, E.K., Hunt, R.J., James, S.R., Foks, N.L., and Komiskey, M.J., 2022, Airborne electromagnetic and magnetic survey data, northeast Wisconsin (ver. 1.1, June 2022): U.S. Geological Survey data release, https://doi.org/10.5066/P93SY9LI.

import matplotlib.pyplot as plt
from os.path import join
from gspy import Survey
from pprint import pprint

Set up the Survey

data_path = '..//..//supplemental//region//WI'
metadata = join(data_path, "data//WI_SkyTEM_survey_md.json")
survey = Survey(metadata)
d_data = join(data_path, 'data//WI_SkyTEM_2021_ContractorData.csv')
d_supp = join(data_path, 'data//WI_SkyTEM_raw_data_md.json')
survey.add_tabular(type='csv', data_filename=d_data, metadata_file=d_supp)

The CRS variable is called spatial_ref and gets initialized in the Survey. The spatial_ref is a dataless coordinate variable, meaning there are no data values, all information is contained within attributes.

print(survey.xarray.spatial_ref)

<xarray.DataArray 'spatial_ref' ()>
array(0.)
Coordinates:
    spatial_ref  float64 0.0
Attributes: (12/18)
    crs_wkt:                           PROJCRS["NAD83(HARN) / Wisconsin Trans...
    semi_major_axis:                   6378137.0
    semi_minor_axis:                   6356752.314140356
    inverse_flattening:                298.257222101
    reference_ellipsoid_name:          GRS 1980
    longitude_of_prime_meridian:       0.0
    ...                                ...
    longitude_of_central_meridian:     -90.0
    false_easting:                     520000.0
    false_northing:                    -4480000.0
    scale_factor_at_central_meridian:  0.9996
    authority:                         EPSG
    wkid:                              3071

The Survey also has a spatial_ref property which returns the spatial_ref variable

print(survey.spatial_ref)

<xarray.DataArray 'spatial_ref' ()>
array(0.)
Coordinates:
    spatial_ref  float64 0.0
Attributes: (12/18)
    crs_wkt:                           PROJCRS["NAD83(HARN) / Wisconsin Trans...
    semi_major_axis:                   6378137.0
    semi_minor_axis:                   6356752.314140356
    inverse_flattening:                298.257222101
    reference_ellipsoid_name:          GRS 1980
    longitude_of_prime_meridian:       0.0
    ...                                ...
    longitude_of_central_meridian:     -90.0
    false_easting:                     520000.0
    false_northing:                    -4480000.0
    scale_factor_at_central_meridian:  0.9996
    authority:                         EPSG
    wkid:                              3071

Grid Mapping

Following the CF conventions on Grid Mappings, the spatial_ref variable should contain key information defining the coordinate reference system. The attribute grid_mapping_name is required. Other key attributes include wkid and crs_wkt.

print('grid_mapping_name: '+survey.xarray.spatial_ref.attrs['grid_mapping_name'])
print('wkid: '+survey.xarray.spatial_ref.attrs['wkid'])
print('crs_wkt: '+survey.xarray.spatial_ref.attrs['crs_wkt'])

grid_mapping_name: transverse_mercator
wkid: 3071
crs_wkt: PROJCRS["NAD83(HARN) / Wisconsin Transverse Mercator",BASEGEOGCRS["NAD83(HARN)",DATUM["NAD83 (High Accuracy Reference Network)",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4152]],CONVERSION["Wisconsin Transverse Mercator 83",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-90,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",520000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-4480000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["State-wide spatial data management."],AREA["United States (USA) - Wisconsin."],BBOX[42.48,-92.89,47.31,-86.25]],ID["EPSG",3071]]

Then, each data variable should have an attribute grid_mapping that references the spatial_ref coordinate variable

pprint(survey.tabular['DEM'].attrs)

{'axis': 'Z',
 'datum': 'North American Vertical Datum of 1988 (NAVD88)',
 'grid_mapping': 'spatial_ref',
 'long_name': 'Digital elevation model',
 'null_value': 'not_defined',
 'positive': 'up',
 'standard_name': 'dem',
 'units': 'meter',
 'valid_range': array([172.52357302, 354.1462611 ])}

Making a new Spatial Ref

If you need to make a new spatial_ref variable, this can be done with GSPy’s Spatial_ref class

from gspy.src.classes.survey.Spatial_ref import Spatial_ref

The Spatial_ref class takes a dictionary of values and looks for a wkid, crs_wkt, or a proj_string in that order. Note, a wkid should have an authority key passed with it either as a separate authority field, or as a colon separated string, e.g., ‘EPSG:4326’. If none is provided EPSG will be used by default.

new_crs = Spatial_ref.from_dict({'wkid': 4326, 'authority': 'EPSG'})

pprint(new_crs)

<xarray.DataArray ()>
array(0.)
Attributes:
    crs_wkt:                      GEOGCRS["WGS 84",ENSEMBLE["World Geodetic S...
    semi_major_axis:              6378137.0
    semi_minor_axis:              6356752.314245179
    inverse_flattening:           298.257223563
    reference_ellipsoid_name:     WGS 84
    longitude_of_prime_meridian:  0.0
    prime_meridian_name:          Greenwich
    geographic_crs_name:          WGS 84
    horizontal_datum_name:        World Geodetic System 1984 ensemble
    grid_mapping_name:            latitude_longitude
    authority:                    EPSG
    wkid:                         4326

Note: If you are resetting the CRS variable in a Survey, be sure that all data groups are also updated to match and all coordinate variables (particularly x, y, and z) need to be updated. In other words, if you change from a projected coordinate system with easting and northing coordinates to a geographic coordinate system, then the x and y coordinate variables need to be changed to longitude and latitude.