Coordinate System

In previous examples the coordinates specified did not correspond to any coordinates in the real world. [0, 0, 0] was an arbitrarily chosen point to represent the origin and the X, Y and Z axes were similarly arbitrarily defined. However, for real world applications it is often useful to define the coordinates based on coordinates in the real world. Coordinate systems allow for the coordinates within the application to map to or from coordinates in the real world.

The Python SDK uses the pyproj library to handle the bulk of its coordinate systems definition. In particular the pyproj.CRS (Coordinate Reference System) class is used to define the common part of the coordinate system. In addition to this, a local transform can be specified to provide slight alterations to the pyproj coordinate system to provide a more convenient coordinate system (for example, to counteract UTM distortion away from the prime meridian).

Note that the Maptek Python SDK only allows for objects to be tagged with a coordinate system. It does not support coordinate system conversions. Adding, removing or changing the coordinate system property of an object has no effect on the topology of the object. If the coordinate system of an object is defined, the Convert Coordinate System menu option in PointStudio (and other Maptek applications) can correctly convert the object to a different coordinate system.

The Maptek Python SDK allows for coordinate systems to be specified in two ways:

The coordinate system property of a Topology object can be assigned to a pyproj.CRS object. The pyproj.CRS object is automatically wrapped in a CoordinateSystem object with no local transform information.

from pyproj import CRS
 
from mapteksdk.project import Project
from mapteksdk.data import Polygon
 
project = Project()
 
with project.new("cad/rectangle_with_coordinate_system", Polygon) as new_polygon:
  # Coordinates should be in latitude, longitude for wgs84.
  new_polygon.points = [[112, -9, 0], [112, -44, 0],
                      [154, -44, 0], [154, -9, 0]]
  new_polygon.coordinate_system = CRS.from_epsg(4326)

Alternatively, the coordinate system of an object can be specified using a CoordinateSystem object. This allows for an additional local transform to be specified via a LocalTransform object. This allows for adjustments to the coordinate system which make it more convenient to use.

import math
 
from pyproj import CRS
 
from mapteksdk.project import Project
from mapteksdk.data import Polygon, CoordinateSystem, LocalTransform
 
project = Project()
 
local_transform = LocalTransform(
  horizontal_origin=[15362089, -4028802],
  horizontal_scale_factor=0.9,
  horizontal_rotation=math.radians(45))
# EPSG 3857 is web mercator. Avoid using this coordinate system
# in real scripts.
coordinate_system = CoordinateSystem(CRS.from_epsg(3857), local_transform)
 
with project.new("cad/rectangle_with_local_transform", Polygon) as new_polygon:
  new_polygon.points = [[12467782, -5465442], [12467782, -1006021],
                      [17143201, -1006021], [17143201, -5465442]]
  new_polygon.coordinate_system = coordinate_system

Note: Note: The coordinate system will not change the points property. This means if there is already a coordinate system set it will not transform the points from the old system to new.

Transforming between two coordinate systems

Maptek recommends using the pyproj library for performing a transform between two coordinate systems. It is recommended that you consult the pyproj documentation for the specifics. However the following is intended as a starting point and quick reference.

The key part is creating a Transformer object from pyproj which has the source coordinate reference system and the target coordinate reference system. From there, call the transform() class method on the transform and provide it with the X and Y coordinates to convert. It returns the transformed coordinates which then need to be written back to the object.

The following is not suitable if a LocalTransform is applied.

from mapteksdk.data import Polygon
from mapteksdk.project import Project
from pyproj import CRS, Transformer
 
project = Project()
with project.new('/cad/Tasmania_WGS84', Polygon) as new_polygon:
    # This is a very loose approximation of the Australian State of mainland
    # of Tasmania.
    new_polygon.points = [
        (144.70488006, -40.50805977, 0),
        (146.41937944, -41.16362087, 0),
        (148.12853241, -40.73132871, 0),
        (148.32398191, -42.23678208, 0),
        (147.83977638, -42.87653518, 0),
        (147.26300257, -43.06258608, 0),
        (146.85094046, -43.63700682, 0),
        (146.03269512, -43.55630696, 0),
        (145.18856871, -42.24195051, 0),
        (145.21375825, -41.98466186, 0),
        (144.60700003, -41.00818821, 0),
        (144.70474642, -40.50794073, 0),
     ]
    new_polygon.coordinate_system = CRS.from_epsg(4326)  # WGS84
 
# Tasmania lies in MGA zone 55.
copied_polygon = project.copy_object(new_polygon.id,
                                     '/cad/Tasmania_MGA_Zone_55')
with project.edit(copied_polygon) as polygon:
    new_coordinate_system = CRS.from_epsg(7855)
    previous_coordinate_system = polygon.coordinate_system.crs
 
    transformer = Transformer.from_crs(previous_coordinate_system,
                                       new_coordinate_system)
 
    polygon.points[:,0], polygon.points[:,1] = transformer.transform(
        polygon.points[:, 1], polygon.points[:, 0])
 
    # The points now refer to the new coordinate system so we update the
    # information stored.
    polygon.coordinate_system = new_coordinate_system

As you see when defining the points, X is the longitude and Y is the latitude however typically when dealing with coordinates, they are normally in the order of latitude and longitude. As a result the above code has to swap the X and Y. See Why is the axis ordering in PROJ not consistent? In PROJ’s Frequently Asked Questions for more.

The above script will produce two polygons as seen below, the change in apparent shape is due to the projection.