Coordinate transformations

Coordinate transformations#

Until now, when we handled multiple shapefiles within one process (e.g., during selections by locations / spatial filters we have mostly worked under the assumption, that both shapefiles are in the same spatial projection. However, in reality this is rarely the case. Instead, we commonly face the situation of having two or more different projections. When we work with QGIS, we often don’t pay attention to this, as there this is handled in the background. When we script/program, we need to pay attention to this ourselves. The first intuitive reaction is then often: “Ok, I use QGIS before to reproject all layers into one common coordinate system, and then I don’t have to worry about this.” This is true, but think about the example we looked at in the context of spatial filters. Would you have taken the WDPA shapefile of >5GB size and reproject it first into another coordinate system? Leaving aside that it probably takes a lot of time to do so, it also produces a lot of additional data that for many cases we are not even interested in (e.g., marine protected areas), or that outside of our area of interest (i.e., the nine countries). So, it makes from many perspectives a lot of sense to think about coordinate transformations “on-the-fly”, and in this chapter we give you some tools on hand to do so. We will work with the same shapefiles as in the spatial filter chapter - with one difference: the gadm-layer is in an equal area projection (i.e., EPSG:3035), and hence in a different coordinate system than the WDPA shapefile:

from osgeo import ogr, osr
folder = "D:/OneDrive - Conservation Biogeography Lab/_TEACHING/__Classes-Modules_HUB/M8_Geoprocessing-with-python/data/"
ds1 = ogr.Open(folder + "EU_vector/gadm36_epsg3035.shp")
ctr = ds1.GetLayer()
ds2 = ogr.Open(folder + "EU_vector/WDPA_sub.shp")
pa = ds2.GetLayer()

If we take now the geometry of the first feature (Germany) and apply the spatial filter, like we did previously, we get:

# Get the geometry
germanyFeat = ctr[0]
germanyGEOM = germanyFeat.geometry()
# Apply the spatial filter
pa.SetSpatialFilter(germanyGEOM)
pa.GetFeatureCount()

0

This differs markedly from what we found previously. In fact, the correct output should be 23196 protected areas intersecting with the borders of Germany. Why is that? Well, let’s have a look at the spatial references of the two layers:

ctr.GetSpatialRef().ExportToWkt()

'PROJCS["ETRS89-extended / LAEA Europe",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","3035"]]'

pa.GetSpatialRef().ExportToWkt()

'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'

Clearly, we see these differences, particularly with regards to the unit of the coordinates, which are meters in case of the country borders, and degrees in case of the protected areas. The task that is ahead of us now is to homogenize these two projections. We can theoretically imagine two ways:

  1. We transform the geometries of all protected areas into EPSG:3035

  2. We transform the geometry of Germany into EPSG:4326 Obviously, transforming one geometry is computationally less intensive that transforming 41158 geometries (i.e., the total number of protected areas inside the shapefile).

What is technically needed to do this, are two things:

  1. Create a ruleset that translates coordinates from one coordinate system to another,

  2. Apply this ruleset to a geometry.

# Create the rule for transforming coordinates from one spatial reference to another
fromCS = ctr.GetSpatialRef()
toCS = pa.GetSpatialRef()
ct = osr.CoordinateTransformation(fromCS, toCS)

Technically, it was not necessary to define the two variables fromCS and toCS. I did it anyways, because this way you remember how the rule is built, and in addition it already looks ahead to other sessions/chapters, e.g., when we do this with a raster and shapefile. In that case, we won’t be able to call the function .GetSpatialRef() for the raster file, as this is a method that is made for vector layers.

# Apply the rule
germanyGEOM.Transform(ct)

0

Now, the geometry of Germany is in the coordinate system of the protected area layer (EPSG:4326), and we can run our spatial filter again:

pa.SetSpatialFilter(germanyGEOM)
pa.GetFeatureCount()

23196

This corresponds now to the expected result. Pretty cool, huh?? And all that without having produced a new shapefile, that occupies disc space and which we won’t need afterwards anymore, as (for this task) we were only interested in the number of protected areas per country. In addition, because coordinate transformations are after all trigonometric calculations, they are very fast. A couple of things you should keep in mind, though:

  1. Doing a coordinate transformation actually manipulates the geometry. This means that you don’t want to work with a reference to the original object (e.g., via .GetGeometryRef()) but to a copy which we obtain via .geometry(). If you want to be sure that you don’t mess up thhe original dataset, you can also create a clone of a geometry. In our case you would do that via germanyGEOM_cl = germanyGEOM.Clone()

  2. You need to think in advance in which direction you want to do the coordinate transformation. The computer will apply any rule to the geometry that you as a user define - but if the rule is wrong, then the outcome will be wrong as well.

Question

Can you think of an operation, where you first transform a geometry into another coordinate system and afterwards back to the original one? How would you make sure here that you don’t mess up the original data file?