Understanding Coordinate Systems

As a geoscience platform, Earth coordinate system support for all spatial data is important. So what are coordinate systems?

Spatial data viewed on a map, or a screen, or represented in a computer database, almost always uses a Cartesian reference frame for coordinate systems. A Cartesian reference frame represents a location by an (x, y, z) coordinate, which is a distance along an X, Y, Z axis in a right-handed sense, with (X,Y) defining the orthogonal horizontal axis, and Z the vertical axis, positive up. An equivalent terminology common with Earth coordinate systems would be (Easting, Northing, Elevation), which makes the understanding of (x, y, z) explicit.

But the Earth is not flat, and for thousands of years cartographers have been constructing maps that illustrate spatial information based on a mapping (a projection) of an actual location on the somewhat elliptical, and bumpy, Earth onto the horizontal plane (X, Y) of a map (and now computer screen). In order for us to know how to different sets of spatial data are located we need to know the details of the Coordinate System of the spatial locations. Many software systems do not deal with this complexity and leave it to the users to ensure that (x, y, z) locations of all data in the project are in the same Coordinate System. Geosoft software, and for that matter all advanced GIS systems, provide extensive support coordinate systems, which as a minimum will provide the ability to fully describe the coordinate system for any set of spatial data.

A complete Coordinate System description has three parts, a horizontal Datum, a Map Projection and a Vertical Datum:

Horizontal Datum

All Coordinate Systems require the definition of the Horizontal Datum (often simply refered to as the Datum), which is the mathematical description of an ellipsoid-like surface that models the actual surface of the Earth. The longitude and latitude of a location is relative to one of the defined Horizontal Datums. Longitude and latitude use a "Geographic Coordinate System", which excludes elevation and not a Cartesian Coordinate System.

Note that the longitude and latitude of a fixed location on earth is relative to a Horizontal Datum. The same fixed location will have a different longitude and latitude for each Horizontal Datum.

There are two basic types Horizontal Datums - geocentric and non-geocentric. WGS 84 has emerged as the global standard geocentric datum, which effectively uses the GRS 80 reference ellipsoid. All modern Coordinate Systems are based either on WGS 84, or a close version that adheres to small nuances for a particular mapping agency. For example, NAD83 (North America) and GDA94 (Australia) are, for the purposes of practical geoscience mapping applications, effectively align with WGS 84 and ses the detic and Geoid. Geocentric datums are locations on a perfect spheroid centered at the geographic center of the Earth, and such Datums have become the standard reference Horizontal Datums equivalent. Older datasets and coordinate systems (generally pre-1980's) use various non–geocentric Horizontal Datums that cartographers established to best fit the area of a particular mapping region. These Horizontal Datums will have a different reference ellipsoid, different earth center, and often a slightly different orientation of the axis. Examples include NAD27 (North America), ED50 (Europe), AGD66 (Australia).

All coordinate systems must have a defined Horizontal Datum. If only the Horizontal Datum is defined coordinates are assumed to be geographic Longitude, Latitude on the defined Horizontal Datum.

Map Projection

Map projections define the mathematics and parameters that transform Longitude Latitude to/from a map Cartesian (x, y) location. There are many different models for doing this, such as Mercator (vertical cylindrical),Transverse Mercator (horizontal cylindrical), Lambert (conic). Every map and spatial dataset that used Cartesian Coordinates will have a defined map projection and parameters. For example, "UTM zone 15N" is a Transverse Mercator projection, with "Universal" implying the standard set of Transverse Mercator parameters for zone 15 North:

origin latitude: 0
central meridian: 93 W
scale factor: 0.9996
false Easting: 500,000
false Northing: 0
distance unit: metre

Vertical Datum

The Vertical Datum is the reference surface from which elevations (the z coordinate) is defined. This surface is either the geodetic datum surface, which is the surface of the Horizontal Datum reference ellipsoid, or it is the geoid, which is gravitational potential surface that most closely matches mean sea level globally. As with Horizontal Datums, modern Vertical Datums are usually based on the geodetic datum surface, while historic elevations are more commonly referred to as "above mean sea level", which means they reference the Earth's geoid.

It is common for Coordinate Systems to not define the Vertical Datum, but for 3D geoscience applications that include an elevation (z), it is important to be explicit about Vertical Datum as the difference between a geodetic elevation and the geoid elevation can be more than 100 m, and is commonly 10's of metres.

In Geosoft, we currently maintain the name of the Vertical Datum as part of the coordinate system descriptive name.

How to define a coordinate system

The simplest way to describe a coordinate system is by using a coordinate system string, which in Geosoft as the following form:

"horizontal_datum / map_projection [vertical_datum]"

Horizontal datum names recognized by Geosoft are listed in reference table: C:\Program Files\Geosoft\Desktop Applications 9\csv\datum.csv. The horizontal_datum, map_projection and vertical_datum should use the well-known names defined by the EPSG Geodetic Parameter Registry.

Map projections recognized by Geosoft are listed in reference table: C:\Program Files\Geosoft\Desktop Applications 9\csv\transform.csv

The vertical_datum string can currently be any descriptive string, including "geoid" or "geodetic", though to future-proof your code we recommend that you use the common short name for a known vertical datum. For example, the "North American Vertical Datum of 1988" is commonly called "NAVD88".

Following are some valid descriptive name examples:

"NAD83 / UTM zone 15N"	horizontal datum: NAD83 map projection: UTM zone 15N vertical datum: not defined
"WGS 84"	horizontal datum: WGS 84 map projection: none, this is a geographic (longitude, latitude) system vertical datum: not defined
"GDA94 / MGA zone 55"	horizontal datum: GDA94 map projection: MGA zone 55, which is equivalent to "UTM zone 55S" vertical datum: not defined
NAD83 [NAVD88]"	horizontal datum: NAD83 map projection: none, this is a geographic (longitude, latitude) system vertical datum: North America Vertical Datum of 1988 (NAVD88)
"NAD83 / UTM zone 15N [geoid]"	horizontal datum: NAD83 map projection: UTM zone 15N vertical datum: 'geoid'

To create an instance of a coordinate system you can provide the coordinate system string to geosoft.gxpy.coordinate_system.Coordinate_system to create a Coordinate_system instance. For example, the following code defines the "UTM zone 15N" projection on the "NAD83" datum, and the "NAD27" geographic coordinate system. These defined systems are then used to translate spatial coordinates from one system to the other.

translate_coordinates_between_systems.py

import geosoft.gxpy.gx as gx
import geosoft.gxpy.coordinate_system as gxcs
import numpy as np

# create context
gxc = gx.GXpy()

# define coordinate systems and a transformer
cs_utm = gxcs.Coordinate_system('NAD83 / UTM zone 15N')
cs_nad27 = gxcs.Coordinate_system('NAD27')
cs_transform = gxcs.Coordinate_translate(cs_utm, cs_nad27)

# example transform a single (x, y) coordinate
lon_lat = cs_transform.convert((345000, 64250000))
print('(lon, lat): {}'.format(lon_lat))

# example transform a single (x, y, elevation) coordinate
print('(lon, lat, elevation): {}'.format(cs_transform.convert((345000, 64250000, 50))))

# example translate a list of (x, y, z) tuples
locations = [(345000, 64250000, 50), (345500, 64250000, 60), (346000, 64250000, 70)]
nad27_locations = cs_transform.convert(locations)
for xyz in nad27_locations:
    print(xyz)

# example transform a numpy array in-place
data = np.array([[345000, 64250000, 50, 55000],
                 [345500, 64250000, 60, 55150],
                 [346000, 64250000, 70, 56000]],
                dtype=float)
nad27_locations = cs_transform.convert(data, in_place=True)
for xyz in data:
    print(xyz)

# compare coordinate systems
print(cs_utm == cs_nad27)
print(gxcs.Coordinate_system('WGS 84') == gxcs.Coordinate_system('WGS 84'))
print(gxcs.Coordinate_system('GDA94 [geodetic]') == gxcs.Coordinate_system('GDA94 [geoid]'))

# output:
# (lon, lat): (88.777242210445195, -38.498998514257273)
# (lon, lat, elevation): (88.777242211469414, -38.498998481053022, 50.0)
# [ 88.77724221 -38.49899848  50.        ]
# [ 88.77151111 -38.49908532  60.        ]
# [ 88.76577999 -38.49917188  70.        ]
# [  8.87772422e+01  -3.84989985e+01   5.00000000e+01   5.50000000e+04]
# [  8.87715111e+01  -3.84990853e+01   6.00000000e+01   5.51500000e+04]
# [  8.87657800e+01  -3.84991719e+01   7.00000000e+01   5.60000000e+04]
# False
# True
# False

line 9,10	Create coordinate systems.
line 11	Create a coordinate system transformer, in this case transforming from "NAD83 / UTM zone 15N" to "NAD27" geographic coordinates.
line 14	Transform an (x, y) location
line 18	Transform a 3D location, (x, y, elevation)
line 22	Transform a list of locations
line 31	Transform a numpy array of locations in-place
line 36	Comparing coordinate systems

Coordinate systems and spatial data

All spatial information in the Geosoft environment will have a defined coordinate system, which if not explicitly set or inherited will be named "*unknown". Geosoft functions that mix spatial data from different coordinate systems will automatically transform coordinates to be in the coordinate system required by the context. Spatial locations that have an "*unknown" coordinate system are generally assumed to match the context in which the coordinates are used.

When data is imported into Geosoft it is good practice to define the coordinate system my assigning the known coordinate system to the coordinate_system property of the data.. If you import data via one of the Geosoft import functions, and the data has a way of describing its coordinate system, the Geosoft import function will make a best effort to set the the coordinate_system property from the data. For example, ESRI data files often have a well defined coordinate system which is recognized when the data is imported.

The geosoft.gxpy module exposes various spatial data structures via classes, such as the Geosoft_database class, the Grid class, the Geometry class, and the View class. All spatial classes in the Geosoft environment will have a property named coordinate_system, which is an instance of the geosoft.gxpy.coordinate_system.Coordinate_system class. The coordinate_system property can be used get or set the instance coordinate systems using any form supported by the Coordinate_system class constructor. The following script demonstrated the versatility of how to establish a coordinate systems using various forms. These examples are shown applied to a geosoft.gxpy.grid.Grid instance, but may be equally applied to any spatial class instance.

assign_coordinate_system.py

import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrd

gxc = gx.GXpy()

# create a memory grid as an example of a spatial object
grid = gxgrd.Grid.new(properties=({'nx': 10, 'ny': 10}))
print(grid.coordinate_system) # initially '*unknown"

# define by a Geosoft-style coordinate system name. Parameters are derived from internal Geosoft tables.
grid.coordinate_system = "NAD83 / UTM zone 17N"
print(grid.coordinate_system)
print(grid.coordinate_system.gxf)
print(grid.coordinate_system.coordinate_dict())

# example use of GXF strings to change the datum to NAD27. Here we remove the name and local datum transform
# and allow the Coordinate_system class to complete parameters for NAD27 from the tables.
gxf = grid.coordinate_system.gxf
gxf[0] = ''
gxf[1] = "NAD27"
gxf[4] = ''
grid.coordinate_system = gxf
print('gxf:', grid.coordinate_system.gxf)

# fully explicit definition of UTM zone 17N on NAD27 datum using GXF strings.
grid.coordinate_system = ['',
                          'NAD27',
                          '"Transverse Mercator",0,-87,0.9996,500000,0',
                          'm,1',
                          '"*local_datum",-8,160,176,0,0,0,0']
print('gxf:', grid.coordinate_system.gxf)

# ... from a json string. Note how to properly escape the string embedded in a string.
js = '{"units": "m,1", "datum": "NAD27", "projection": "\\"Transverse Mercator\\",0,-87,0.9996,500000,0"}'
grid.coordinate_system = js
print('json:', grid.coordinate_system.gxf)

# ... from an ESRI WKT string
wkt = 'PROJCS["NAD_1927_UTM_Zone_16N",' + \
          'GEOGCS["GCS_North_American_1927",' + \
          'DATUM["D_North_American_1927",' + \
          'SPHEROID["Clarke_1866",6378206.4,294.9786982]],' + \
          'PRIMEM["Greenwich",0.0],' + \
          'UNIT["Degree",0.0174532925199433]],' + \
          'PROJECTION["Transverse_Mercator"],' + \
          'PARAMETER["False_Easting",500000.0],' + \
          'PARAMETER["False_Northing",0.0],' + \
          'PARAMETER["Central_Meridian",-87.0],' + \
          'PARAMETER["Scale_Factor",0.9996],' + \
          'PARAMETER["Latitude_Of_Origin",0.0],' + \
          'UNIT["Meter",1.0],' + \
          'AUTHORITY["EPSG",26716]]'
grid.coordinate_system = wkt
print('from wkt:', grid.coordinate_system.esri_wkt)

Running this script will print the following:

*unknown
NAD83 / UTM zone 17N
['NAD83 / UTM zone 17N', 'NAD83,6378137,0.0818191910428158,0', '"Transverse Mercator",0,-81,0.9996,500000,0', 'm,1', '"NAD83 to WGS 84 (1)",0,0,0,0,0,0,0']
{'orientation': '', 'type': 'Geosoft', 'vcs': '', 'units': 'm,1', 'name': 'NAD83 / UTM zone 17N', 'local_datum': '"NAD83 to WGS 84 (1)",0,0,0,0,0,0,0', 'datum': 'NAD83,6378137,0.0818191910428158,0', 'projection': '"Transverse Mercator",0,-81,0.9996,500000,0'}
gxf: ['NAD27 / UTM zone 17N', 'NAD27,6378206.4,0.0822718542230039,0', '"Transverse Mercator",0,-81,0.9996,500000,0', 'm,1', '"NAD27 to WGS 84 (4)",-8,160,176,0,0,0,0']
gxf: ['NAD27 / UTM zone 16N', 'NAD27,6378206.4,0.0822718542230039,0', '"Transverse Mercator",0,-87,0.9996,500000,0', 'm,1', '*local_datum,-8,160,176,0,0,0,0']
json: ['NAD27 / UTM zone 16N', 'NAD27,6378206.4,0.0822718542230039,0', '"Transverse Mercator",0,-87,0.9996,500000,0', 'm,1', '"NAD27 to WGS 84 (4)",-8,160,176,0,0,0,0']
from wkt: PROJCS["NAD_1927_UTM_Zone_16N",GEOGCS["GCS_North_American_1927",DATUM["D_North_American_1927",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-87.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",26716]]