Python – Rotating Projected Coordinates Points Horizontally

affine transformationcoordinate systempython

I have this Sentinel-2 picture that uses the WGS84 UTM zone 23S as projection type. My picture is rescaled to 3413 by 4741. My objective is to transform the X and Y pixels to EPSG 32723 and then from that to 4326.
I'm trying to achieve that by doing an affine transformation with the Cartesian product of my width range and height range. My affine transformation uses the following data: 20metres(that equals to 1 pixel), 0 degrees, -20metres, 441519.5898437500 for X and 7988872.0410156259
for Y. Code below:

pixel_transform = affine.Affine(
        20,
        0,
        441509.58984375,       
        0,
        -20, 
       7988882.041015626
)

After that, before converting the X and Y points to degrees longitude and latitude, I create a GeoPandas dataframe to make a shape file and check in QGIS if my coordinates are getting transformed properly.
My issue is: apparently everything is being properly transformed as my coordinates points are being projected in QGIS in the exact shape of the real world location they represent EXCEPT for the fact that everything is kinda of rotating anti-clockwise by somewhat 90 degrees and I can't really grasp why. I also tried to validate my coordinates using Kepler.GL (this time converting from 32723 to 4326) and it also rotates in the same way. So my points, although they "draw" the correct shape of what I'm projecting, don't fit in the map, as they for some reason rotate.
I assume I'm missing something during my affine transformation, but I really don't know what. I'm using a pre-processed Sentinel-2 photo, so maybe some important data got lost in the way and I should add some degrees to the transformation? I really don't know. Has anyone experienced something similar?

Best Answer

I figured it out. Turns out my image, probably due to the preprocessing, had it's X and Y axis inverted. Then to get my proper coordinates, I needed to switch the elements in my affine matrix.

Related Solutions

[GIS] Transforming data from a rotated pole lat-lon grid into regular lat-lon coordinates in QGIS

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

    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:

[GIS] Convert Lat/Lon sexagesimal degrees to UTM coordinates

The conversion you did was correct, as in, your lat/lon in ETRS89 30N does indeed translate to (1301616.242, 4173536.779).

The problem is in the zone you are using. UTM zones go from -3º to +3º from your central meridian, which is set at 500km, in the X axis. At 37º of latitude, each 1º is about 88.5km, meaning the entire zone spans from 234.5km to 765.5km. If your point is beyond that, it's time to change zones. Which is exactly what is happening.

Your derived coordinates at zone 30N have an X of 1301616.242m, that is way more than what's acceptable, and you should change zones. Converting to zone 31N, you get an X of 768365m, which is pretty close to the value that this other web conversor used. Hence, the conversor probably identified the correct zone, and translated accordingly.

One other thing that is to note is that, as user Marco mentioned in the comments, it may be that your longitude is wrong, as in, it should be -6º (6º west) where it's written 6º. This way, using Zone 30N would be correct (incidently, your X in this setup would be about 230525m).

Related Question