Spatial Projections

(This is based entirely on the GDAL/OSR Tutorial and Python GDAL/OGR Cookbook.)

The ArchGDAL.SpatialRef, and ArchGDAL.CoordTransform types are lightweight wrappers around GDAL objects that represent coordinate systems (projections and datums) and provide services to transform between them. These services are loosely modeled on the OpenGIS Coordinate Transformations specification, and use the same Well Known Text format for describing coordinate systems.

Coordinate Systems

There are two primary kinds of coordinate systems. The first is geographic (positions are measured in long/lat) and the second is projected (such as UTM - positions are measured in meters or feet).

  • Geographic Coordinate Systems: A Geographic coordinate system contains information on the datum (which implies an spheroid described by a semi-major axis, and inverse flattening), prime meridian (normally Greenwich), and an angular units type which is normally degrees.

  • Projected Coordinate Systems: A projected coordinate system (such as UTM, Lambert Conformal Conic, etc) requires and underlying geographic coordinate system as well as a definition for the projection transform used to translate between linear positions (in meters or feet) and angular long/lat positions.

Creating Spatial References

spatialref = ArchGDAL.importEPSG(2927)
Spatial Reference System: +proj=lcc +lat_0=45.3333333333333 + ... _defs
ArchGDAL.toPROJ4(spatialref)
"+proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5 +lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000.0001016 +y_0=0 +ellps=GRS80 +units=us-ft +no_defs"

The details of how to interpret the results can be found in http://proj4.org/usage/projections.html.

In the above example, we constructed a SpatialRef object from the EPSG Code 2927. There are a variety of other formats from which SpatialRefs can be constructed, such as

We currently support a few export formats too:

Reprojecting a Geometry

source = ArchGDAL.importEPSG(2927)
Spatial Reference System: +proj=lcc +lat_0=45.3333333333333 + ... _defs
target = ArchGDAL.importEPSG(4326)
Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs
ArchGDAL.createcoordtrans(source, target) do transform
    point = ArchGDAL.fromWKT("POINT (1120351.57 741921.42)")
    println("Before: $(ArchGDAL.toWKT(point))")
    ArchGDAL.transform!(point, transform)
    println("After: $(ArchGDAL.toWKT(point))")
end
Before: POINT (1120351.57 741921.42)
After: POINT (47.3488070138318 -122.598149943144)

Reprojecting from a layer

# Plotting with native GEOJSON geographic CRS
p_WGS_84 = AG.getfeature(layer, 0) do feature
    AG.getgeom(feature, 0) do geom
        plot(geom; fa=0.1, title="WGS 84")
    end
end

# Plotting with local projected CRS
p_Lambert_93 = AG.getfeature(layer, 0) do feature
    AG.getgeom(feature, 0) do geom
        source = AG.getspatialref(geom)
        target = AG.importEPSG(2154)
        AG.createcoordtrans(source, target) do transform
            plot(AG.transform!(geom, transform); fa=0.1, title="Lambert 93")
        end
    end
end

plot(p_WGS_84, p_Lambert_93; size=(600, 200), layout=(1,2))