geodeticengine 1.0.14

Library for transformation and conversion of coordinates used in Equinor.

Geodetic Engine

A Python library using the Equinor Geodetic Engine API and pyproj to transform coordinates between different Coordinate Reference Systems (CRS).

Installation

1. Access to the Equinor Geodetic Engine API is required, which can be obtained through AccessIT.

2. A personal subscription key is also necessary to authenticate through APIM. Sign in to api.equinor.com and subscribe to Enterprise.

3. As usual the package it´s installed running pip install:

pip install geodeticengine

Authentication

There are two ways to authenticate:

• User access - MSAL Device code flow
• Application access - MSAL Client credential flow

User access - MSAL Device code flow

This package uses the PublicClientApplication Python class from the MSAL library for user authentication. Bearer tokens are retrieved using the acquire_token_interactive() method, which can be accessed via a local browser on devices that support a Graphical User Interface (GUI). If a device lacks GUI support, such as GitHub Code Spaces, the initiate_device_flow() method generates a specific URL for the user to visit and follows a standard authentication and login process.

In order to use the package with user access to production, you will only need to add one environment variable to your system:

EGE_SUBS_KEYS=<your-subscription-key-for-geodeticengine-prod>

In order to use the package with user access to test, you will need to add the following environment variables to your system:

EGE_API_ENV=test
EGE_SUBS_KEYS=<your-subscription-key-for-geodeticengine-test>

EGE_SUBS_KEYS: This variable holds your APIM subscription key for each environment (prod and test). This will allow the package to access the resources you have permission to use.
EGE_API_ENV: This variable is used to access the Geodetic Engine Test environment. If this environment variable is not set, the package will use the production environment by default. It can also be set to production, by EGE_API_ENV=prod.

User access - Token cache

The token for each environment is cached and stored in the user's home directory, eliminating the need to authenticate before every session. Although an access token expires after one hour, the presence of a cached Refresh Token allows a new Access Token to be obtained without requiring re-authentication. The Refresh Token lasts for 90 days, then you have to log in again.

Application access - MSAL Client credential flow

In order to use the package with application access to production, you will need to add the following environment variables to your system:

EGE_CLIENT_IDS=<your-app-id-prod>
EGE_CLIENT_SECRETS=<your-app-secret-prod>
EGE_SUBS_KEYS=<your-subscription-key-for-geodeticengine-prod>

In order to use the package with application access to test:

EGE_API_ENV=test
EGE_CLIENT_IDS=<your-app-id-prod>;<your-app-id-test>
EGE_CLIENT_SECRETS=<your-app-secret-prod>;<your-app-secret-test>
EGE_SUBS_KEYS=<your-subscription-key-for-geodeticengine-prod>;<your-subscription-key-for-geodeticengine-test>

EGE_API_ENV: This variable is used to access the different environments test and production. If this environment variable is not set, the package will use the production environment by default. It can also be set to prod, by EGE_API_ENV=prod. EGE_CLIENT_IDS: This variable holds your application (client) ID for each environment.
EGE_CLIENT_SECRETS: This variable is used for your application's client secret for each environment. If this variable is not set, the package will automatically fall back to user access for authentication.
EGE_SUBS_KEYS: This variable holds your APIM subscription key for each environment (prod and test). This will allow the package to access the resources that you have permission to use.

Note that if the EGE_CLIENT_IDS, EGE_CLIENT_SECRETS and EGE_SUBS_KEYS environment variables hold the IDs and secrets for both the production and test environments, the order is important. The production values must always come before the test values.

Transformation grids

Transformation grids are essential for achieving high accuracy when performing datum transformations. The package will automatically fetch required grid if it doesn't already exist on your local machine.

More information on the data available is located in pyproj documentation.

Transformation and Conversion

Transform list of lists :

Transforming and converting coordinates provided as a list of lists can be performed in two ways using this Python package. The first method involves passing all the required information, including the coordinates to be transformed, to the API. In this case, the transformation is carried out in the cloud by the API itself. The second approach is to use the API to request the transformation pipeline—a string containing all the parameters needed to transform or convert the coordinates from CRS A to CRS B. Then, you can use this pipeline to transform the coordinates locally using the transform_coordinates_with_pipeline function. This latter approach is more efficient and well-suited for large datasets, as it does not require sending the coordinates over the network.

Transform geometries:

To transform geometries, there are two possible approaches. You can create a CoordTrans object and then use the method transform_geometry to carry out the transformation, or you can use the standalone function transform_geometry_with_pipeline if you already have a transformation pipeline. The transformation will happen "locally" in both approaches. The geometries will not be sent to the API.

Code Examples

1. Transformation examples

   • Transforming a geometry defined by WKT using the instance of the CoordTrans class
   • Transforming a geometry defined by WKB (with base64 decoding) using the instance of the CoordTrans class
   • Transforming geometry defined by WKB as a byte string with a predefined transformation pipeline
   • Transforming geometry defined by WKT and returning the transformed geometry as WKB represented by a byte string
   • Convert from WGS84 UTM Zone 31N to WGS84 Geographic 2D (EPSG:4326)
   • Transforming coordinates from an engineering CRS to a global CRS and back to engineering again
   • Transform coordinates by specifying EPSG codes for both the source and target CRS and the transformation between them and by sending the coordinates to the API
   • Transform coordinates using a bound CRS and by sending the coordinates to the API
   • Transform coordinates without sending them to the API. Only the transformation pipeline is extracted from the API and the transformation is done locally using pyproj
   • Transform coordinates by using two different transformations
   • Convert geographic coordinates to projected coordinates
   • Transform coordinates by using numpy array as input, and return the transformed coordinates as a numpy array

2. CRS and CTs listing

   • Get all bound CRS to WGS84
   • Get all bound Projected and Projected CRS
   • Get all transformations between ED50 and WGS 84
   • List all available transformations between ED50 and WGS84

Transformation examples

Transforming a geometry defined by WKT using the instance of the CoordTrans class

from geodeticengine import CoordTrans

geom_wkt = "POLYGON((4.4545121 59.553535, 4.555555 60.66666666, 3.5844187 61.11111111, 4.8888888 58.15441217, 4.4545121 59.553535))"
crs_from = "(EPSG:4230, EPSG:1613)"
crs_to = "ST_WGS84_UTM31N_P32631"

transformed_geom = CoordTrans(crs_from=crs_from, crs_to=crs_to, geometry=geom_wkt).transform_geometry(output_fmt="wkt")

Transforming a geometry defined by WKB with base64 decoding using the instance of the CoordTrans class

from geodeticengine import CoordTrans

geom_wkb = "AQYAAAACAAAAAQMAAAABAAAABAAAAEWrSab5YRFAjT1ikTLwTUAmfYx/wNMRQJkT7nRr/k1Ab3zU/IRFEkCNPWKRMvBNQEWrSab5YRFAjT1ikTLwTUABAwAAAAEAAAAEAAAAt0FQ4jG3FkCAZ9at+eFNQAHD0dnyJxdAjT1ikTLwTUBbnIyFfI8XQIBn1q354U1At0FQ4jG3FkCAZ9at+eFNQA=="
crs_from = "ST_ED50_S62_T1613"
crs_to = "ST_WGS84_UTM31N_P32631"

transformed_geom = CoordTrans(crs_from=crs_from, crs_to=crs_to, geometry=geom_wkb).transform_geometry(output_fmt="wkt", num_decimals=3)

Transforming geometry defined by WKB as a byte string with a predefined transformation pipeline

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_geometry_with_pipeline

geom_wkb = b'\x01\x06\x00\x00\x00\x02\x00\x00\x00\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\x9a\xa9\xe9\x17\x1e\x8e!A]\xc4$\xd0\xc7RYAw\x127+\x95\xbc!A\xc4>)\xdb\xfd^YA\x04\x87 bL\xef!A>\xa2M;\x0cSYA\x9a\xa9\xe9\x17\x1e\x8e!A]\xc4$\xd0\xc7RYA\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\x00\xcc\xea\x03\x00\xd9#A\xbd\xd6\xc0\xef\xeaHYAN\xf1\xe5;:\x05$AE<jy>UYA\x010\xa9\xd8\xa45$A\xa5\x81\x92\x86gIYA\x00\xcc\xea\x03\x00\xd9#A\xbd\xd6\xc0\xef\xeaHYA'
crs_from = "ST_WGS84_UTM31N_P32631"
crs_to = "ST_ED50_S62_T1613"

# Use the instanse of the CoordTrans class to retrieve transformation pipeline from Geodetic Engine
ct_pipeline = CoordTrans(crs_from=crs_from, crs_to=crs_to).get_pipeline()

# Use the pipeline requested from the Geodetic Engine to transform the geometry "locally"
transformed_geom = transform_geometry_with_pipeline(geom_wkb, ct_pipeline, output_fmt="wkt", num_decimals=8)

Transforming geometry defined by WKT and returning the transformed geometry as WKB represented by a byte string

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_geometry_with_pipeline

geom_wkb = b'\x01\x06\x00\x00\x00\x02\x00\x00\x00\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\x9a\xa9\xe9\x17\x1e\x8e!A]\xc4$\xd0\xc7RYAw\x127+\x95\xbc!A\xc4>)\xdb\xfd^YA\x04\x87 bL\xef!A>\xa2M;\x0cSYA\x9a\xa9\xe9\x17\x1e\x8e!A]\xc4$\xd0\xc7RYA\x01\x03\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\x00\xcc\xea\x03\x00\xd9#A\xbd\xd6\xc0\xef\xeaHYAN\xf1\xe5;:\x05$AE<jy>UYA\x010\xa9\xd8\xa45$A\xa5\x81\x92\x86gIYA\x00\xcc\xea\x03\x00\xd9#A\xbd\xd6\xc0\xef\xeaHYA'
crs_from = "ST_WGS84_UTM31N_P32631"
crs_to = "ST_ED50_S62_T1613"

# Use the instanse of the CoordTrans class to retrieve transformation pipeline from Geodetic Engine
ct_pipeline = CoordTrans(crs_from=crs_from, crs_to=crs_to).get_pipeline()

# Use the pipeline requested from the Geodetic Engine to transform the geometry "locally"
transformed_geom = transform_geometry_with_pipeline(geom_wkb, ct_pipeline, output_fmt="wkt", num_decimals=8)

Convert from WGS84 UTM Zone 31N to WGS84 Geographic 2D

from geodeticengine import CoordTrans
from geodeticengine.geospatial_utils import plot_geometry_on_map

geom_wkt = geom_wkt = 'POLYGON ((443793.60417810833 6629500.057140259, 449427.01243148284 6997077.285665373, 595811.6738373992 6998103.368741313, 606520.5686072728 6630605.5272681555, 443793.60417810833 6629500.057140259))'
crs_from = "ST_WGS84_UTM31N_P32631"
crs_to = "EPSG:4326"

# Convert from WGS84 UTM Zone 31N to WGS84 Geographic 2D (EPSG:4326)
transformed_geom = CoordTrans(crs_from=crs_from, crs_to=crs_to, geometry=geom_wkt).transform_geometry(return_object=True)

Transforming coordinates from an engineering CRS to a global CRS and back to engineering again

from geodeticengine import CoordTrans

points = [[8911.832,5139.165]]
crs_from = "EQUINOR:4100002" # ST_MON_Refinery_E4100002
crs_to = "EPSG:25832"        # ST_ETRS89_UTM32N_P25832
ct_from = "EQUINOR:3100002"  # ST_MON_ETRS89_P25832_T3100002
ct = CoordTrans(crs_from=crs_from, crs_to=crs_to, ct_from=ct_from, points=points)
transformed_coords = ct.transform_pointlist()
print(f"Coordinates transformed to {crs_to}:\n{transformed_coords}")

# Transform back to engineering
ct = CoordTrans(crs_from=crs_to, crs_to=crs_from, ct_to=ct_from, points=transformed_coords)
transformed_coords_back = ct.transform_pointlist()
print(f"Coordinates transformed back to the engineering CRS {crs_to}:\n{transformed_coords_back}")

Transform coordinates by specifying EPSG codes for both the source and target CRS and the transformation between them and by sending the coordinates to the API

from geodeticengine import CoordTrans

points = [[8911.832,5139.165]]
crs_from = "EQUINOR:4100002" # ST_MON_Refinery_E4100002
crs_to = "EPSG:25832"        # ST_ETRS89_UTM32N_P25832
ct_from = "EQUINOR:3100002"  # ST_MON_ETRS89_P25832_T3100002
ct = CoordTrans(crs_from=crs_from, crs_to=crs_to, ct_from=ct_from, points=points)
transformed_coords = ct.transform_pointlist()
print(f"Coordinates transformed to {crs_to}:\n{transformed_coords}")

# Transform back to engineering
ct = CoordTrans(crs_from=crs_to, crs_to=crs_from, ct_to=ct_from, points=transformed_coords)
transformed_coords_back = ct.transform_pointlist()
print(f"Coordinates transformed back to the engineering CRS {crs_to}:\n{transformed_coords_back}")

Transform coordinates using a bound CRS and by sending the coordinates to the API

from geodeticengine import CoordTrans

points = [[9,65],[12,70]]
crs_from = "ST_ED50_T1133"
crs_to = "ST_WGS84_G4326"

# Transform coordinates
ct = CoordTrans(crs_from=crs_from, crs_to=crs_to, points=points)
print(f"Transformed coordinates:{ct.transform_pointlist()}")

Transform coordinates without sending them to the API. Only the transformation pipeline is extracted from the API and the transformation is done locally using pyproj

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_coordinates_with_pipeline

crs_from = "EPSG:32632" # WGS84 / UTM zone 32N
crs_to = "EPSG:23031"   # ED50 / UTM zone 31N
ct_code = "EPSG:1613"   # Transformation from  ED50 til WGS84

points = [[273178.571, 6846980.063]]

# Create a transformer object of the source crs, target crs and the transformation between them and extract pipeline
trans_pipeline = CoordTrans(crs_from, crs_to, ct_from=ct_code).get_pipeline()

# Use the transformation pipeline extracted in the line above to transform the coordiantes
transformed_coord_lst = transform_coordinates_with_pipeline(points, trans_pipeline)
print(f"Transformed coordinates by using the pipeline from the GeodeticEngine:\n{transformed_coord_lst}")

Transform coordiantes by using two different transformations

First transform the coordinates from ETRS89 to WGS84 using the EQUINOR:3000034 transformation, and then transform from WGS84 to ED50 using EPSG:1613

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_coordinates_with_pipeline

crs_from = "EPSG:25832"          # ETS89 / UTM zone 32N
crs_to = "EPSG:23031"            # ED50 / UTM zone 31N
ct_from_code = "EQUINOR:3000034" # Equinor custom transformation from ETRS89 to WGS84
ct_to_code = "EPSG:1613"         # Transformation from ED50 to WGS84

points = [[273178.571, 6846980.063]]

# Create a transformer object of the source crs, target crs and the transformation between them
trans_obj = CoordTrans(crs_from, crs_to, ct_from=ct_from_code, ct_to=ct_to_code)

# Get the transformation pipeline
trans_pipeline = trans_obj.get_pipeline(pretty_fmt=True)
print(f"The transformation pipeline for the two transformations:\n{trans_pipeline}\n")


# Use the transformation pipeline extracted in the line above to transform the coordiantes
transformed_coord_lst = transform_coordinates_with_pipeline(points, trans_pipeline)
print(f"Transformed coordinates by using the pipeline from the GeodeticEngine:\n{transformed_coord_lst}")

Convert geographic coordinates to projected coordinates

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_coordinates_with_pipeline

crs_from = "EPSG:4326"           # WGS84 geograpich 2D (lat/lon)
crs_to = "EPSG:32631"            # WGS84 / UTM zone 31N

points = [[1.97926020, 60.53835841]] # longitude first, then latitude (degree)

# Create a transformer object of the source crs, target crs
trans_obj = CoordTrans(crs_from, crs_to)

# Get the transformation pipeline
conv_pipeline = trans_obj.get_pipeline()
print(f"The pipeline for the conversion:\n{trans_pipeline}\n")

# Use the transformation pipeline extracted in the line above to convert the coordinates
converted_coord_lst = transform_coordinates_with_pipeline(points, conv_pipeline)
print(f"Converted coordinates by using the pipeline from the GeodeticEngine:\n{converted_coord_lst}")

Tranform coordinates by using numpy array as input, and return the transformed coordinates as a numpy array

from geodeticengine import CoordTrans
from geodeticengine.CoordTrans import transform_coordinates_with_pipeline
import numpy as np

crs_from = "EPSG:32632" # WGS84 / UTM zone 32N
crs_to = "EPSG:23031"   # ED50 / UTM zone 31N
ct_code = "EPSG:1613"   # Transformation from  ED50 til WGS84

points = np.array([[273178.571, 6846980.063]])

# Create a transformer object of the source crs, target crs and the transformation between them and extract pipeline
trans_pipeline = CoordTrans(crs_from, crs_to, ct_from=ct_code).get_pipeline()

# Use the transformation pipeline extracted in the line above to transform the coordiantes
transformed_coord = transform_coordinates_with_pipeline(points, trans_pipeline, return_ndarray=True)
print(f"Transformed coordinates by using the pipeline from the GeodeticEngine:\n{transformed_coord}")

CRS and CTs listing examples

Get all bound CRS to WGS84

from geodeticengine import CrsSearch
import json

# Get all bound CRS to WGS84
crs_query = CrsSearch(types=["bound projected","bound geographic 2D", "bound dynamic geographic 2D"], target_crs="ST_WGS84_G4326")
all_bound_crs = json.dumps(crs_query.get_entities(), indent=2)
print(f"All bound CRS:\n{all_bound_crs}")

Get all bound Projeced and Projected CRS

from geodeticengine import CrsSearch
import json

crs_query = CrsSearch(types=["bound projected","projected"])
all_bound_proj_and_proj = json.dumps(crs_query.get_entities(), indent=2)
print(f'All "Bound Projeced" and "Projected":\n{all_bound_proj_and_proj}')

Get all transformations between ED50 and WGS 84

from geodeticengine import CtSearch
import json

ct_query = CtSearch(types=["transformation","concatenated operation"], source_crs="ST_ED50_G4230", target_crs="ST_WGS84_G4326")
cts = json.dumps(ct_query.get_entities(), indent=2)
print(f'All transformation between ED50 and WGS 84:\n{cts}')

List all available transformations between ED50 and WGS84

from geodeticengine import CtSearch
import json

ct_query = CtSearch(types=["transformation","concatenated operation"], source_crs="ST_ED50_G4230", target_crs="ST_WGS84_G4326")
all_available_trans = json.dumps(ct_query.get_entities(), indent=2)
print(f'All availiable transformation between ED50 and WGS 84 in total:\n{all_available_trans}')