[GIS] Finding correct extent of projection

coordinate systemextentsprojreprojection-mathematics

It is tricky to manage the extent of a coordinate system when that extent is not an aligned rectangle in another coordinate system.

For example consider Thursday Island, located at:

  • 142.2°, -10.2° (WGS84, EPSG:4326) or equivalently,
  • 1146km, -1156.5km (Australian Albers projection, EPSG:3577).

According to SpatialReference.org (similar values to epsg.io) this latter projection has:

  • WGS84 Bounds: 108.0000, -45.0000, 155.0000, -10.0000
  • Projected Bounds: -2690013.3995, -5098960.4467, 2579169.8548, -1281018.4757

Notice that Thursday Island is well inside of the lat/lon bounds, however it (along with sections of Australian mainland at both Cape York Peninsula and Kakadu) falls outside of the projected bounds.

Say I want to create a web map application using OpenLayers, I need to know the maximum extent of the map in the projected coordinates. The problem is that if I source this information from a standard online reference, I will mistakenly cut off parts of Australia even when choosing a projection that was specifically designed to cover Australia.

Presumably the root cause is that the extent has been represented by an aligned bounding-box and then, rather than transforming the entire boundary of this region precisely (before generating a new bounding box), only two vertices of the box have been properly projected (the bottom-left and top-right corners). Such a method is invalid due to curvature/nonlinearity.

How is this problem normally handled? Is there any authorative source for the extent of a map projection? In particular, is Proj.4 able to report the maximum domain of validity for a transformation? I'd like a solution that works in all cases, including those where the region of validity does not resemble a rectangular coordinate patch (e.g. rHEALPix and QSC).

Best Answer

The primary source of these extents is the EPSG Geodetic Parameter Registry, e.g.:

EPSG Geodetic Parameter Registry

PROJ version 6+

If you have a newer version of PROJ (version 6 or later), use projinfo to extract a WKT2_2018 string:

$ projinfo -o WKT2_2018 EPSG:3577
WKT2_2018 string:
PROJCRS["GDA94 / Australian Albers",
    BASEGEOGCRS["GDA94",
        DATUM["Geocentric Datum of Australia 1994",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4283]],
    CONVERSION["Australian Albers",
        METHOD["Albers Equal Area",
            ID["EPSG",9822]],
        PARAMETER["Latitude of false origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",132,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",-18,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",-36,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Australia - onshore"],
        BBOX[-43.7,112.85,-9.86,153.69]],
    ID["EPSG",3577]]

See the USAGE at the bottom for AREA and BBOX.

pyproj

Using pyproj >=2.0, you can access various properties of a CRS in Python:

In [1]: from pyproj import CRS                                                                                                                      

In [2]: crs = CRS.from_epsg(3577)                                                                                                                   

In [3]: crs                                                                                                                                         
Out[3]: 
<Projected CRS: EPSG:3577>
Name: GDA94 / Australian Albers
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Australia - onshore
- bounds: (112.85, -43.7, 153.69, -9.86)
Coordinate Operation:
- name: Australian Albers
- method: Albers Equal Area
Datum: Geocentric Datum of Australia 1994
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

In [4]: crs.area_of_use                                                                                                                             
Out[4]: AreaOfUse(name=Australia - onshore, west=112.85, south=-43.7, east=153.69, north=-9.86)

In [5]: crs.area_of_use.bounds                                                                                                                      
Out[5]: (112.85, -43.7, 153.69, -9.86)

EPSG.io

Websites like EPSG.io (code) have converted the EPSG registry to a SQLite db (gml.sqlite), which is queried to extract the XML and show the same information. You can also fetch a "bbox" from their API, e.g. http://epsg.io/?q=3577&format=json

"bbox": [-9.86, 112.85, -43.7, 153.69]

Note that these bounds are different than SpatilaReference.org, shown in the question, but that's because that website was based on an older edition of the EPSG database. I generally advise to not use SpatilaReference.org as it is essentially abandoned.

Polygon extents

If bounding boxes are not enough, you can register a free account and download "Area polygons in Shapefile format", which has many overlapping areas:

world

zooming into the point of interest, highlighting urn:ogc:def:area:EPSG::2575 - Australia - onshore:

australia

This is 1 of nearly 20 different areas defined for different coordinate reference systems at that location.

Related Question