Coordinate Systems

In scenarios where real-world coordinates are not necessary, an arbitrary coordinate system can be used, where [0, 0, 0] represents the origin, and the X, Y, and Z axes are defined arbitrarily. However, for applications that require real-world accuracy, it is advantageous to use a defined coordinate system. Coordinate systems enable mapping within the application to real-world coordinates, ensuring consistency and accuracy in spatial data representation.

The Maptek Python SDK employs the pyproj library to handle most of its coordinate system definitions. Specifically, the pyproj.CRS (Coordinate Reference System) class is used to define the core aspects of the coordinate system. Additionally, a local transform can be applied to adjust the pyproj coordinate system slightly, making it more practical for specific use cases, such as correcting UTM distortion away from the prime meridian.

Important:  The Maptek Python SDK only permits tagging objects with a coordinate system. It does not support coordinate system conversions directly. Changing the coordinate system property of an object does not alter the object's topology. If an object’s coordinate system is defined, the Convert Coordinate System tool in PointStudio and other Maptek applications can convert the object to a different coordinate system.

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

  • Direct assignment: 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)
    
    
  • With local transform: 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 method allows for modifications to the coordinate system to suit specific needs.

    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, 0], [12467782, -1006021, 0],
                          [17143201, -1006021, 0], [17143201, -5465442, 0]]
      new_polygon.coordinate_system = coordinate_system
    
    

    Note:  Assigning a coordinate system does not modify the object’s points. This means if there is already a coordinate system set, it will not transform the points from the previous system to the new one.

Transforming between two coordinate systems

We recommend using the pyproj library to transform coordinates between two coordinate systems. Consult the pyproj documentation for detailed guidance. The following steps provide a basic reference:

  1. Create a Transformer object from pyproj, specifying both the source and target coordinate reference systems.

  2. Use the transform() method of the Transformer object to convert the X and Y coordinates.

The key part is creating a Transformer object from pyproj that 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 script demonstrates transforming from one coordinate system to another:

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

Note:  The approach used in this script is not suitable if a local transform is applied.

When defining coordinates, X typically represents longitude and Y represents latitude. However, coordinates are often given in the order of latitude and longitude. Consequently, you may need to swap the X and Y values in the code. For more information on axis ordering in PROJ, see Why is the axis ordering in PROJ not consistent? in the PROJ Frequently Asked Questions.

The script above generates two polygons, as shown below. The change in apparent shape is due to the projection applied.