# Transforming data from a rotated pole lat-lon grid into regular lat-lon coordinates in QGIS?

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.

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.

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. Jul 3 '17 at 21:07

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

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")

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])

#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:

• The GRASS7 algorithm `v.mkgrid` does basically the same. 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. 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. Jul 4 '17 at 19:32

If you want to convert this for CORDEX data so you can use a free online converter: Rotation of Coordinates Based On CORDEX Domains. But if your database is not CORDEX data so you can use Open NC File for any netcdf data. In other hands, maybe this link would be useful: Can Open NC File to open a NARCCAP NetCDF files

• please be sure to mention your interest in the solution you are promoting Dec 28 '19 at 18:12