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 SpatialRef
s can be constructed, such as
ArchGDAL.importEPSG(::Int)
: based on the EPSG codeArchGDAL.importEPSGA(::Int)
: based on the EPSGA codeArchGDAL.importESRI(::String)
: based on ESRI projection codesArchGDAL.importPROJ4(::String)
: based on the PROJ.4 string (reference)ArchGDAL.importURL(::String)
: download from a given URL and feed it intoSetFromUserInput
for you.ArchGDAL.importWKT(::String)
: WKT stringArchGDAL.importXML(::String)
: XML format (GML only currently)
We currently support a few export formats too:
ArchGDAL.toMICoordSys(spref)
: Mapinfo style CoordSys format.ArchGDAL.toPROJ4(spref)
: coordinate system in PROJ.4 format.ArchGDAL.toWKT(spref)
: nicely formatted WKT string for display to a person.ArchGDAL.toXML(spref)
: converts into XML format to the extent possible.
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)
using DataFrames
import GeoFormatTypes as GFT
coords = zip(rand(10), rand(10))
df = DataFrame(geom=ArchGDAL.createpoint.(coords), name="test");
df.geom = ArchGDAL.reproject(df.geom, GFT.EPSG(4326), GFT.EPSG(28992))
10-element Vector{ArchGDAL.IGeometry{ArchGDAL.wkbPoint}}: Geometry: POINT (-579409.221301706 -5619443.91672546) Geometry: POINT (-582336.53377174 -5641651.45511224) Geometry: POINT (-570196.229261112 -5721836.35931894) Geometry: POINT (-511734.935279046 -5692122.02899001) Geometry: POINT (-472709.692881904 -5730761.23840706) Geometry: POINT (-526881.251439192 -5700632.80878514) Geometry: POINT (-448858.365704017 -5629142.66755675) Geometry: POINT (-462901.900742049 -5705682.18556205) Geometry: POINT (-502511.080998322 -5733997.21354294) Geometry: POINT (-538791.633578617 -5722901.51503972)
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))
References
Some background on OpenGIS coordinate systems and services can be found in the Simple Features for COM, and Spatial Reference Systems Abstract Model documents available from the Open Geospatial Consortium. The GeoTIFF Projections Transform List may also be of assistance in understanding formulations of projections in WKT. The EPSG Geodesy web page is also a useful resource. You may also consult the OGC WKT Coordinate System Issues page.