Tutorial for the Geo Add-on

This tutorial shows how to load a GPS track into a geo located DXF file and also the inverse operation, exporting geo located DXF entities as GeoJSON files.

Please read the section Intended Usage in the documentation of the ezdxf.addons.geo module first.

Warning

TO ALL BEGINNERS!

If you are just learning to work with geospatial data, using DXF files is not the way to go! DXF is not the first choice for storing data for spatial data analysts. If you run into problems I cannot help you as I am just learning myself.

The complete source code and test data for this tutorial are available in the github repository:

https://github.com/mozman/ezdxf/tree/master/docs/source/tutorials/src/geo

Setup Geo Location Reference

The first step is setting up the geo location reference, which is not doable with ezdxf yet - this feature may come in the future - but for now you have to use a CAD application to do this. If the DXF file has no geo location reference the projected 2D coordinates are most likely far away from the WCS origin (0, 0), use the CAD command “ZOOM EXTENDS” to find the data.

Load GPX Data

The GPX format stores GPS data in a XML format, use the ElementTree class to load the data:

def load_gpx_track(p: Path) -> Iterable[Tuple[float, float]]:
    """Load all track points from all track segments at once."""
    gpx = ET.parse(p)
    root = gpx.getroot()
    for track_point in root.findall(".//gpx:trkpt", GPX_NS):
        data = track_point.attrib
        # Elevation is not supported by the geo add-on.
        yield float(data["lon"]), float(data["lat"])

The loaded GPS data has a WSG84 EPSG:4326 projection as longitude and latitude in decimal degrees. The next step is to create a GeoProxy object from this data, the GeoProxy.parse() method accepts a __geo_interface__ mapping or a Python object with a __geo_interface__ attribute/property. In this case as simple “LineString” object for all GPS points is sufficient:

def add_gpx_track(msp, track_data, layer: str):
    geo_mapping = {
        "type": "LineString",
        "coordinates": track_data,
    }
    geo_track = geo.GeoProxy.parse(geo_mapping)

Transform the data from the polar representation EPSG:4326 into a 2D cartesian map representation EPSG:3395 called “World Mercator”, this is the only projection supported by the add-on, without the need to write a custom transformation function:

    geo_track.globe_to_map()

The data is now transformed into 2D cartesian coordinates in meters and most likely far away from origin (0, 0), the data stored in the GEODATA entity helps to transform the data into the DXF WCS in modelspace units, if the DXF file has no geo location reference you have to stick with the large coordinates:

    # Load geo data information from the DXF file:
    geo_data = msp.get_geodata()
    if geo_data:
        # Get the transformation matrix and epsg code:
        m, epsg = geo_data.get_crs_transformation()
    else:
        # Identity matrix for DXF files without a geo location reference:
        m = Matrix44()
        epsg = 3395
    # Check for compatible projection:
    if epsg == 3395:
        # Transform CRS coordinates into DXF WCS:
        geo_track.crs_to_wcs(m)
        # Create DXF entities (LWPOLYLINE)
        for entity in geo_track.to_dxf_entities(dxfattribs={"layer": layer}):
            # Add entity to the modelspace:
            msp.add_entity(entity)
    else:
        print(f"Incompatible CRS EPSG:{epsg}")

We are ready to save the final DXF file:

doc.saveas(str(out_path))

In BricsCAD the result looks like this, the underlying images were added by the BricsCAD command MAPCONNECT and such a feature is not planned for the add-on:

../_images/gpx_tracks.png

Export DXF Entities as GeoJSON

This will only work with a proper geo location reference, the code shown accepts also WCS data from DXF files without a GEODATA object, but the result is just unusable - but in valid GeoJSON notation.

First get epsg code and the CRS transformation matrix:

    # Get the geo location information from the DXF file:
    geo_data = msp.get_geodata()
    if geo_data:
        # Get transformation matrix and epsg code:
        m, epsg = geo_data.get_crs_transformation()
    else:
        # Identity matrix for DXF files without geo reference data:
        m = Matrix44()

Query the DXF entities to export:

    for track in msp.query("LWPOLYLINE"):
        export_geojson(track, m)

Create a GeoProxy object from the DXF entity:

def export_geojson(entity, m):
    # Convert DXF entity into a GeoProxy object:
    geo_proxy = geo.proxy(entity)

Transform DXF WCS coordinates in modelspace units into the CRS coordinate system by the transformation matrix m:

    # Transform DXF WCS coordinates into CRS coordinates:
    geo_proxy.wcs_to_crs(m)

The next step assumes a EPSG:3395 projection, everything else needs a custom transformation function:

    # Transform 2D map projection EPSG:3395 into globe (polar)
    # representation EPSG:4326
    geo_proxy.map_to_globe()

Use the json module from the Python standard library to write the GeoJSON data, provided by the GeoProxy.__geo_interface__ property:

    # Export GeoJSON data:
    name = entity.dxf.layer + ".geojson"
    with open(TRACK_DATA / name, "wt", encoding="utf8") as fp:
        json.dump(geo_proxy.__geo_interface__, fp, indent=2)

The content of the GeoJSON file looks like this:

{
  "type": "LineString",
  "coordinates": [
    [
      15.430999,
      47.06503
    ],
    [
      15.431039,
      47.064797
    ],
    [
      15.431206,
      47.064582
    ],
    [
      15.431283,
      47.064342
    ],
    ...
}

Custom Transformation Function

This sections shows how to use the GDAL package to write a custom transformation function. The example reimplements the builtin transformation from unprojected WGS84 coordinates to 2D map coordinates EPSG:3395 “World Mercator”:

from osgeo import osr
from ezdxf.math import Vec3

# GPS track in WGS84, load_gpx_track() code see above
gpx_points = list(load_gpx_track('track1.gpx'))

# Create source coordinate system:
src_datum = osr.SpatialReference()
src_datum.SetWellKnownGeoCS('WGS84')

# Create target coordinate system:
target_datum = osr.SpatialReference()
target_datum.SetWellKnownGeoCS('EPSG:3395')

# Create transformation object:
ct = osr.CoordinateTransform(src_datum, target_datum)

# Create GeoProxy() object:
geo_proxy = GeoProxy.parse({
    'type': 'LineString',
    'coordinates': gpx_points
})

# Apply a custom transformation function to all coordinates:
geo_proxy.apply(lambda v: Vec3(ct.TransformPoint(v.x, v.y)))

The same example with the pyproj package:

from pyproj import Transformer
from ezdxf.math import Vec3

# GPS track in WGS84, load_gpx_track() code see above
gpx_points = list(load_gpx_track('track1.gpx'))

# Create transformation object:
ct = Transformer.from_crs('EPSG:4326', 'EPSG:3395')

# Create GeoProxy() object:
geo_proxy = GeoProxy.parse({
    'type': 'LineString',
    'coordinates': gpx_points
})

# Apply a custom transformation function to all coordinates:
geo_proxy.apply(lambda v: Vec3(ct.transform(v.x, v.y)))

Polygon Validation by Shapely

Ezdxf tries to avoid to create invalid polygons from HATCH entities like a hole in another hole, but not all problems are detected by ezdxf, especially overlapping polygons. For a reliable and robust result use the Shapely package to check for valid polygons:

import ezdxf
from ezdxf.addons import geo
from shapely.geometry import shape

# Load DXF document including HATCH entities.
doc = ezdxf.readfile('hatch.dxf')
msp = doc.modelspace()

# Test a single entity
# Get the first DXF hatch entity:
hatch_entity = msp.query('HATCH').first

# Create GeoProxy() object:
hatch_proxy = geo.proxy(hatch_entity)

# Shapely supports the __geo_interface__
shapely_polygon = shape(hatch_proxy)

if shapely_polygon.is_valid:
    ...
else:
    print(f'Invalid Polygon from {str(hatch_entity)}.')

# Remove invalid entities by a filter function
def validate(geo_proxy: geo.GeoProxy) -> bool:
    # Multi-entities are divided into single entities:
    # e.g. MultiPolygon is verified as multiple single Polygon entities.
    if geo_proxy.geotype == 'Polygon':
        return shape(geo_proxy).is_valid
    return True

# The gfilter() function let only pass compatible DXF entities
msp_proxy = geo.GeoProxy.from_dxf_entities(geo.gfilter(msp))

# remove all mappings for which validate() returns False
msp_proxy.filter(validate)

Interface to GDAL/OGR

The GDAL/OGR package has no direct support for the __geo_interface__, but has builtin support for the GeoJSON format:

from osgeo import ogr
from ezdxf.addons import geo
from ezdxf.render import random_2d_path
import json

p = geo.GeoProxy({'type': 'LineString', 'coordinates': list(random_2d_path(20))})
# Create a GeoJSON string from the __geo_interface__ object by the json
# module and feed the result into ogr:
line_string = ogr.CreateGeometryFromJson(json.dumps(p.__geo_interface__))

# Parse the GeoJSON string from ogr by the json module and feed the result
# into a GeoProxy() object:
p2 = geo.GeoProxy.parse(json.loads(line_string.ExportToJson()))