Python API Spatial Reference Systems

Creating Spatial Reference Systems

In IlwisObjects you have several different ways to create Coordinate Systems. First you can just load a Coordinate System from an existing csy file, by simply setting the working catalog to the right folder and passing the name of the file as an argument to the constructor of the Coordinate System.

Engine.setWorkingCatalog("file:///D:/<working>/<directory>/")

csy = CoordinateSystem("Cochabamba.csy")

But of course there are also ways to create a new Coordinate System if you don't have a csy file containing all the information yet. One way is by using the epsg number of the according spatial reference system you want to use.

csy = CoordinateSystem("code=epsg:23035")

As you can see the way to pass the epsg number is to make it a string and prepend code=. The coordinate system now has the spatial reference system that is used for Europe - 24°E to 30°E and ED50 by country.

Another way of loading a Spatial Reference System is by using a custom proj4 definition. Through this method you can create a highly customized Spatial Reference System for your personal needs.

csy = CoordinateSystem("code=proj4:+proj=utm +zone=35 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")

The projection now has the name 'utm', is in the utm zone 35, using the intl ellipsoid, a 7 term datum transformation to WGS84 and has meters as unit. For more detailed information on all possible parameters you can have a look at this wiki site of Proj.4.

Transformation

With the help of IlwisObjects it is also possible to transform geometries from one Coordinate System to another. It is simply done by calling the .transform() method on a gemoetry and passing it the Coordinate System to which it should be transformed.

pm = CoordinateSystem("code=epsg:3035")
wgs = CoordinateSystem("code=epsg:4326")
proj4 = CoordinateSystem("code=proj4:+proj=utm +zone=35 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")

g = Geometry("POINT(766489.647 6840642.671)", pm)

g2 = g.transform(wgs)
# POINT (6.8854936501947668 52.2235232557146460)
g3 = g.transform(proj4)
# POINT (-866479.8875676055904478 5979085.6121951946988702)