affine_transform xy coords from gda94

Question

I'm trying to figure out how to convert a polygon whose coordinates are in Spatial Reference GDA94 (EPSG 4283) into xy coordinates (inverse affine transformation matrix).

The following code works:

import sys

import numpy as np

from osgeo import gdal
from osgeo import gdalconst

from shapely.geometry import Polygon
from shapely.geometry.polygon import LinearRing

# Bounding Box (via App) approximating part of QLD.
poly = Polygon(
    LinearRing([
        (137.8, -10.6),
        (153.2, -10.6),
        (153.2, -28.2),
        (137.8, -28.2),
        (137.8, -10.6)
    ])
)

# open raster data
ds = gdal.Open(sys.argv[1], gdalconst.GA_ReadOnly)

# get inverse transform matrix
(success, inv_geomatrix) = gdal.InvGeoTransform(ds.GetGeoTransform())
print inv_geomatrix

# build numpy rotation matrix
rot = np.matrix(([inv_geomatrix[1], inv_geomatrix[2]], [inv_geomatrix[4], inv_geomatrix[5]]))
print rot

# build numpy translation matrix
trans = np.matrix(([inv_geomatrix[0]], [inv_geomatrix[3]]))
print trans

# build affine transformation matrix
affm = np.matrix(([inv_geomatrix[1], inv_geomatrix[2], inv_geomatrix[0]],
                  [inv_geomatrix[4], inv_geomatrix[5], inv_geomatrix[3]],
                  [0, 0, 1]))
print affm

# poly is now a shapely geometry in gd94 coordinates -> convert to pixel
# - project poly onte raster data
xy = (rot * poly.exterior.xy + trans).T  # need to transpose here to have a list of (x,y) pairs

print xy

Here's the output of the printed matrices:

(-2239.4999999999995, 20.0, 0.0, -199.49999999999986, 0.0, -20.0)
[[ 20.   0.]
 [  0. -20.]]
[[-2239.5]
 [ -199.5]]
[[  2.00000000e+01   0.00000000e+00  -2.23950000e+03]
 [  0.00000000e+00  -2.00000000e+01  -1.99500000e+02]
 [  0.00000000e+00   0.00000000e+00   1.00000000e+00]]
[[ 516.5   12.5]
 [ 824.5   12.5]
 [ 824.5  364.5]
 [ 516.5  364.5]
 [ 516.5   12.5]]

Is there a way to do this with scipy.ndimage's affine_transform function?

score 1 · Accepted Answer · edited May 23 '17 at 12:04

There are a few options. Not all spatial transformations are in linear space, so they can't all use an affine transform, so don't always rely on it. If you have two EPSG SRIDs, you can do a generic spatial transform with GDAL's OSR module. I wrote an example a while back, which can be adapted.

Otherwise, an affine transform has basic math:

                    / a  b xoff \ 
[x' y' 1] = [x y 1] | d  e yoff |
                    \ 0  0   1  /
or
    x' = a * x + b * y + xoff
    y' = d * x + e * y + yoff

which can be implemented in Python over a list of points.

# original points
pts = [(137.8, -10.6),
       (153.2, -10.6),
       (153.2, -28.2),
       (137.8, -28.2)]

# Interpret result from gdal.InvGeoTransform
# see http://www.gdal.org/classGDALDataset.html#af9593cc241e7d140f5f3c4798a43a668
xoff, a, b, yoff, d, e = inv_geomatrix

for x, y in pts:
    xp = a * x + b * y + xoff
    yp = d * x + e * y + yoff
    print((xp, yp))

This is the same basic algorithm used in Shapely's shapely.affinity.affine_transform function.

from shapely.geometry import Polygon
from shapely.affinity import affine_transform

poly = Polygon(pts)

# rearrange the coefficients in the order expected by affine_transform
matrix = (a, b, d, e, xoff, yoff)

polyp = affine_transform(poly, matrix)
print(polyp.wkt)

Lastly, it's worth mentioning that the scipy.ndimage.interpolation.affine_transform function is intended for image or raster data, and not vector data.

Yes the shapely.affinity module implements this nicely and I can use it to transform the polygons in EPSG 4283 (GDA94) to xy coordinates for indexing into a numpy (GeoTiff) array (also in GDA94). This doesn't really answer the question of whether this can be done by converting the polygon to a numpy array and performing an affine_transform on it directly (avoiding for loops and pure python which might be slower for polygons that contain > 500,000 points). — James Mills, Apr 16 '13 at 00:19
I think my colleague and I worked out that scipy.ndimage.affine_transform is designed ot work with images and not vectors/polygons. The best way to achieve the general solution is to perform an affine_transform such as your solution above; which shapely.affinity.affine_transform implements. I will accept your answer if you revise it to say "it can't be done with scipy.ndimage" and that you should perform a vector/polygon affine transform (liek yours) or use shapely.affinity.affine_transform --James — James Mills, Apr 16 '13 at 00:29

affine_transform xy coords from gda94

1 Answers1