Spatial Reference

Spatial reference systems define a geodetic date (i.e. defined by an ellipsoid like WGS84 or GRS80) and a projection to define a transformation between 3D coordinates and a flat surface. There are three different encodings of a spatial reference system, which allow an easier exchange of such information and to assign them as properties of geo-referenced data formats like GeoTIFF or NetCDF:

• EPSG: This is the most compact declaration of a spatial reference system. An integer is assigned to a spatial reference system, which is often the preferred solution in many GIS interfaces (https://epsg.org/home.html).

• PROJ4/PROJ6: PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another. It comes along with an API, which translates “PROJ-Strings” (e.g., “+proj=merc +lat_ts=56.5 +ellps=GRS80”) into a representation of coordinate reference system (https://proj.org/).

• WKT: WKT stands for “well-known text representation of coordinate reference systems” and is a low-level, human-readable string format for spatial reference systems, similar to “PROJ-Strings”, but with every information already included within the string. WKT will become deprecated and make way for WKT2, which standardises and defines some new attributes.

geospade provides a class SpatialRef inside its crs module, which is a very slim interface toward all different spatial reference system representations. Below we define three different encodings for the same spatial reference system, namely “MGI GK M34”:

from geospade.crs import SpatialRef

epsg_code = 31259

wkt_string = 'PROJCS[\"MGI / Austria GK M34\",GEOGCS[\"MGI\",DATUM[\"Militar_Geographische_Institute\",SPHEROID' \
             '[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],' \
             'TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],AUTHORITY[\"EPSG\",\"6312\"]],' \
             'PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,' \
             'AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4312\"]],PROJECTION[\"Transverse_Mercator\"],' \
             'PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",16.33333333333333],' \
             'PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",750000],' \
             'PARAMETER[\"false_northing\",-5000000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],' \
             'AUTHORITY[\"EPSG\",\"31259\"]]'

proj4_params = '+init=EPSG:31259 +proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=750000 +y_0=-5000000 +ellps=bessel ' \
               '+towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs'

The constructor of SpatialRef has two arguments: the encoded spatial reference system (EPSG, PROJ4 or WKT) and an optional argument for pre-defining the type (“epsg”, “wkt” or “proj4”). If no type is specified, the type of the spatial reference system is guessed:

# initialise spat. ref. system by specifying the type
sref = SpatialRef(epsg_code, sref_type="epsg")
# automatically guessing the type if no type is provided
sref = SpatialRef(epsg_code)

This initialised object now provides direct access to all three spatial reference system represenations.

sref.epsg

31259

sref.wkt

'PROJCS["MGI / Austria GK M34",GEOGCS["MGI",DATUM["Militar_Geographische_Institute",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],AUTHORITY["EPSG","6312"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4312"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",16.33333333333333],PARAMETER["scale_factor",1],PARAMETER["false_easting",750000],PARAMETER["false_northing",-5000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","31259"]]'

sref.proj4

'+proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=750000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs'

For each attribute, a check is performed if the forward conversion with respect to the spatial reference type used in the constructor is successful. In the case above, a transformation from EPSG to PROJ4 has worked, but not the other way around:

sref_proj4 = SpatialRef(sref.proj4)
sref_proj4.epsg

C:\Users\cnavacch\AppData\Local\Continuum\miniconda3\envs\geospade\lib\site-packages\geospade-0.0.post0.dev66+gf2d2bab-py3.6.egg\geospade\crs.py:521: UserWarning: Conversion from 'PROJ4' to 'EPSG' is not possible.
  warnings.warn(warn_msg.format(self._sref_type.upper(), tar_sref_type.upper()))

Moreover, some additional properties are given to work with commonly used Python objects. This includes an osgeo.osr.SpatialReference object,

sref.osr_sref

<osgeo.osr.SpatialReference; proxy of <Swig Object of type 'OSRSpatialReferenceShadow *' at 0x000001DF7BAFEF90> >

a pretty WKT string,

print(sref.to_pretty_wkt())

PROJCS["MGI / Austria GK M34",
    GEOGCS["MGI",
        DATUM["Militar_Geographische_Institute",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232],
            AUTHORITY["EPSG","6312"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4312"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",16.33333333333333],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",750000],
    PARAMETER["false_northing",-5000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","31259"]]

a PROJ4 parameter dictionary,

sref.to_proj4_dict()

{'proj': 'tmerc',
 'lat_0': 0.0,
 'lon_0': 16.33333333333333,
 'k': 1.0,
 'x_0': 750000.0,
 'y_0': -5000000.0,
 'ellps': 'bessel',
 'towgs84': '577.326,90.129,463.919,5.137,1.474,5.297,2.4232',
 'units': 'm',
 'no_defs': True}

and finally a Cartopy projection system.

sref.to_cartopy_proj()

C:\Users\cnavacch\AppData\Local\Continuum\miniconda3\envs\geospade\lib\site-packages\cartopy\mpl\feature_artist.py:163: UserWarning: Unable to determine extent. Defaulting to global.
  warnings.warn('Unable to determine extent. Defaulting to global.')

<cartopy.crs.TransverseMercator object at 0x000001DF79F1BAF0>

Moreover, you can initialise a SpatialRef object from a osgeo.osr.SpatialReference object.

from osgeo import osr

osr_sref = osr.SpatialReference()
osr_sref.ImportFromEPSG(4326)
sref_osr = SpatialRef.from_osr(osr_sref)
sref_osr.proj4

'+proj=longlat +datum=WGS84 +no_defs'

Several static methods can also help you transform between different representations, e.g.

SpatialRef.osr_to_wkt(osr_sref)

'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]'

More SpatialRef objects can also interact with each on the basis of comparison operators:

sref == sref_osr

False

sref != sref_osr

True