6

I am working with data from a numerical weather model, the output of which is given in grid cells where I have the lat-lon coordinates of each cell's center point. The lat-lon specifications are, however, in a rotated coordinate system with the south pole at (lon: 10, lat: -40). The problem I need help with is how one transforms a grid (a vector layer) which is defined in the rotated lat-lon coordinates into the regular lat-lon coordinates inside QGIS. The figure below shows my specific problem of a rectangular grid in pole rotated lat-lon coordinates (off the coast of West Africa) and the target points in regular lat-lon coordinates over Denmark. Notice that the square grid doesn’t simply get translated (moved) and rotated, there is also a lattitude-dependent scaling issue.

Square grid in pole rotated lat-lon coordinates (at the Equator at West Africa's coast) needs to be transformed to a grid in regular lat-lon (here over Denmark).

A previous thread on SE shows how the transformation can be calculated for points with the help of some simple trigonometry (Lon/lat transformation). These calculations work perfectly (for points) when I implement them in R, but I need to do this transformation for a vector layer inside QGIS. It is common for meteorological offices to only specify the coordinates of the South or North Pole in the rotated lat-lon grid. The QGIS-based transformation therefore has to be done the way simondk does it in his reply to the thread I referred to (he implemented a short piece of Matlab code that does this here.

I've tried to do the transformation with the Affine Transformation plugin where I calculate the parameters for the plugin based on this previous thread. This procedure almost gets it right - the grid is rotated to the right general location (above Denmark), but there seems to be something wrong with the scaling. The center points of the transformed grid cells are not exactly at the correct location - this is clearly seen in the corners of the grid. In the image below, my transformed grid is red, while the correct grid center point locations are the black crosses.

I'm not sure what goes wrong here. Can Affine Transformation even do this lat-lon rotation properly?

I've done extensive Googling and there seem to be many people out there with similar problems, but no solutions. Here is a link to a dropbox folder that contains the data I’ve used for the figures, in case you want to test a solution on the actual rotated grid and target points I’m working with.

Black crosses are the correct locations for the center points, while the grid is my post-transformation result.

I prefer QGIS.

  • If you wish to also ask about ArcGIS Desktop then you can always do that as a separate question. – PolyGeo Jul 3 '17 at 20:23
  • Rotation can be afforded with PyQGIS code. – xunilk Jul 3 '17 at 21:07
4

I created a 0.5 x 0.5 grid (degrees) at your interest area:

enter image description here

and developed next PyQGIS code (by using formulas in your link) for arbitrarily rotate this grid 70º around its center (9.71943416148, 55.4053054573).

from math import sin, cos, pi

def rotatePoints(axis, r_angle, points):

    points_rot = []

    for point in points:
        tr = QgsPoint( point.x() - axis.x(), point.y() - axis.y() )
        xrot = tr.x()*cos(r_angle*(pi/180)) - tr.y()*sin(r_angle*(pi/180))  
        yrot = tr.x()*sin(r_angle*(pi/180)) + tr.y()*cos(r_angle*(pi/180)) 
        xnew = xrot + axis.x()
        ynew = yrot + axis.y()
        points_rot.append(QgsPoint(xnew,ynew))

    return points_rot

def createPolygonMemoryLayer(epsg, pol_rot):

    uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

    mem_layer = QgsVectorLayer(uri,
                               "Polygon",
                               "memory")

    QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

    prov = mem_layer.dataProvider()

    n = len(pol_rot)

    feats = [ QgsFeature() for i in range(n) ]

    for i, feat in enumerate(feats):
        feat.setGeometry(QgsGeometry.fromPolygon(pol_rot[i]))
        feat.setAttributes([i])
        mem_layer.addFeature(feat, True)

    prov.addFeatures(feats)

#Code starts here

layer = iface.activeLayer()

crs = layer.crs()
epsg = crs.postgisSrid()

features = layer.getFeatures()

geom = []

for feature in features:
    geom.append(feature.geometry().asPolygon())

points = geom

n = layer.featureCount()

pol_rot = [ [] for i in range(n) ]

axis = QgsPoint(9.71943416148, 55.4053054573)

r_angle = 70

for i in range(n):
    points_rot = rotatePoints(axis, r_angle, points[i][0])
    pol_rot[i].append(points_rot)

createPolygonMemoryLayer(epsg, pol_rot)

After running above code at the Python Console of QGIS, I got a new rotated grid without any distortion; as it can be observed at next image:

enter image description here

  • The GRASS7 algorithm v.mkgrid does basically the same. – AndreJ Jul 4 '17 at 14:50
  • Thanks xunilk for the quality reply! Based on your answer I have edited my question with some more information. I’m fairly green in Python and haven’t used PyQGIS before but as I understand your code, it is able to rotate a grid around a point. The pole rotation that a lot of meteorological data is provided in has a scaling component to it as well though (in the edits I’ve made to the question, I have specified this further). But I think you are close to home! My guess is that your rotatePoints function only needs a small tweak, I’m just not sure how to implement it myself in Python. – jonaswp Jul 4 '17 at 19:32
  • Rather than taking an axis point and a rotation angle as arguments, it should just take the coordinates of the south pole in the rotated grid (this is the only information that Met Offices provide about their transformation – we never have a simple rotation angle). I’ll share a link to a short bit of Matlab code that does the job here: se.mathworks.com/matlabcentral/fileexchange/… I’ve implemented an R-version of this earlier based on that Matlab function, and it was very straightforward to do. – jonaswp Jul 4 '17 at 19:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.