James Mills
James Mills

Reputation: 19050

affine_transform xy coords from gda94

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?

Upvotes: 3

Views: 1625

Answers (1)

Mike T
Mike T

Reputation: 43682

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.

Upvotes: 1

Related Questions