Geographical frames (`pycraf.geospatial`)#

Introduction#

GIS software often works with coordinate frames other than GPS/WGS84. Therefore, the geospatial sub-package offers routines to convert between different geographical projections/frames. Some of the common use cases would be the transformations between GPS/WGS84 and UTM, ETSR89, or ITRF. However, the geospatial sub-package provides a factory to produce arbitrary transforms, based on EPSG/ESRI codes or proj4 strings.

Note

Internally, the pycraf.geospatial module uses the pyproj package, which itself is a wrapper around the proj.4 software for the transformation. In principle, one could use pyproj directly, but pycraf.geospatial offers unit checking (by making use of the units functionality.)

There are thousands of pre-defined frames (which one can work with by simply using the correct code, e.g., EPSG or ESRI; see Wikipedia), but the underlying proj4 library also accepts user-defined frames. For example, the ETRS89 frame is defined with the following:

+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
+ellps=GRS80 +units=m +no_defs

Using `pycraf.geospatial`#

The provided functions are just a thin wrapper around pyproj, mostly improving convenience for the user:

>>> import pycraf.geospatial as geo
>>> import astropy.units as u

>>> # coordinates of the 100-m telescope at Effelsberg/Germany
>>> rt_lon, rt_lat = 6.88361 * u.deg, 50.52483 * u.deg

>>> etrs_lon, etrs_lat = geo.wgs84_to_etrs89(rt_lon, rt_lat)
>>> print(etrs_lon, etrs_lat)  
4100074.484338885 m 3050584.470522307 m

If the built-in conversions are not sufficient, you can simply build your own transform with the transform_factory. For example, if one wants Gauss-Kruger coordinates (for the Germany tile):

>>> transform = geo.transform_factory(
...     geo.geospatial.EPSG.WGS84, 31467
...     )
>>> transform(rt_lon, rt_lat)  
(<Quantity 3350002.4611806367 m>, <Quantity 5600924.810577951 m>)

The factory even produces a correct doc-string:

>>> print(transform.__doc__)

    Convert coordinates from `sys_in` to `sys_out`.

    - `sys_in`: EPSG.WGS84
    - `sys_out`: EPSG:31467

    Parameters
    ----------
    lon, lat : `~astropy.units.Quantity`
        `sys_in` longitude and latitude [deg]

    Returns
    -------
    x, y : `~astropy.units.Quantity`
        `sys_out` world coordinates [m]

The EPSG and ESRI enums contain some pre-defined projections for convenience.

Note

There are many projections that only produce accurate results over a relatively small region. In some cases such as UTM, there is a bulk of pre-defined frames, one for each so-called zone. Since UTM is quite often used in GIS applications, geospatial has a helper function utm_zone_from_gps to get the correct zone name for a certain Geographical position. Likewise, with epsg_from_utm_zone one could directly query the EPSG code (if needed):

>>> utm_zone = geo.utm_zone_from_gps(rt_lon, rt_lat)
>>> utm_zone
'32N'

>>> utm_lon, utm_lat = geo.wgs84_to_utm(rt_lon, rt_lat, utm_zone)
>>> print(utm_lon, utm_lat)  
349988.58854241937 m 5599125.388981281 m

Imperial units#

The proj4 and pyproj software packages allow to work with frames that are historical or of regional interest, only. Some of these don’t work with units of Meters, but feet, etc. Per default, pyproj (versions prior to 2.0) converts all world (=physical) units to Meters (input and/or output). One can specifically ask for the original units by doing:

>>> import pycraf.geospatial as geo
>>> import astropy.units as u
>>> from astropy.units.imperial import ft
>>> import pyproj

>>> proj_wgs84 = pyproj.Proj('epsg:4326')
>>> # Louisiana South (ftUS)
>>> proj_nad83 = pyproj.Proj('epsg:3452', preserve_units=False)

>>> my_trafo = pyproj.Transformer.from_crs(
...     proj_wgs84.crs, proj_nad83.crs, always_xy=True
...     )
>>> my_trafo.transform(-92.105819, 30.447921)  
(925806.5486332772, 216168.1432314818)

This is the wrong result. But pyproj.Proj has an option:

>>> proj_nad83 = pyproj.Proj('epsg:3452', preserve_units=True)

that gives the correct result:

>>> my_trafo = pyproj.Transformer.from_crs(
...     proj_wgs84.crs, proj_nad83.crs, always_xy=True
...     )
>>> my_trafo.transform(-92.105819, 30.447921)  
(3037416.9849743457, 709211.6499186204)

The geospatial sub-package makes life a bit easier, because we can use the units conversion:

>>> transform = geo.transform_factory(4326, 3452)
>>> x, y = transform(-92.105819 * u.deg, 30.447921 * u.deg)
>>> x, y  
(<Quantity 925806.5486332772 m>, <Quantity 216168.1432314818 m>)

>>> x.to(ft), y.to(ft)  
(<Quantity 3037416.9849743457 ft>, <Quantity 709211.6499186204 ft>)

>>> transform = geo.transform_factory(3452, 4326)
>>> transform(3037416.985 * ft, 709211.650 * ft)  
(<Quantity -92.10581908573734 deg>, <Quantity 30.447921938477027 deg>)

Unfortunately, there seems to be no way to ask pyproj about which original unit is tied to an EPSG code (although it internally must know it, otherwise, it wouldn’t work correctly).

Geocentric systems#

Some frames are Geocentric systems that work with (x, y, z) coordinates. One important example is ITRF:

>>> geo.wgs84_to_itrf2005(rt_lon, rt_lat, 300 * u.m)  
(<Quantity 4033874.035684814 m>,
 <Quantity 486981.61007349996 m>,
 <Quantity 4900340.557846253 m>)

Reference/API#

pycraf.geospatial Package#

This subpackage contains a couple of functions to convert coordinates between typical GIS frames, such as GPS (WGS84), UTM, and ETRS89. The underlying transformations are based on the pyproj package, which itself is a wrapper around the proj.4 software.

Functions#

`epsg_from_utm_zone`(utm_zone)	EPSG number for UTM zone code.
`etrs89_to_wgs84`(elon, elat)	Convert ETSR89 coordinates to GPS/WGS84.
`itrf2005_to_wgs84`(x, y, z)	Convert ITRF2005 coordinates to GPS/WGS84.
`itrf2008_to_wgs84`(x, y, z)	Convert ITRF2008 coordinates to GPS/WGS84.
`transform_factory`(sys_in, sys_out[, ...])	A factory to produce conversion functions (decorated with `ranged_quantity_input` for unit handling).
`utm_to_wgs84`(ulon, ulat, utm_zone)	Convert UTM coordinates to GPS/WGS84.
`utm_zone_from_gps`(glon, glat)	UTM zone code from GPS coordinates.
`wgs84_to_etrs89`(glon, glat)	Convert GPS/WGS84 coordinates to ETSR89.
`wgs84_to_itrf2005`(glon, glat, height)	Convert GPS/WGS84 coordinates to ITRF2005.
`wgs84_to_itrf2008`(glon, glat, height)	Convert GPS/WGS84 coordinates to ITRF2008.
`wgs84_to_utm`(glon, glat, utm_zone)	Convert GPS/WGS84 coordinates to UTM.

Classes#

`EPSG`(value[, names, module, qualname, type, ...])	Enum with often used EPSG codes.
`ESRI`(value[, names, module, qualname, type, ...])	Enum with often used ESRI codes.