Reputation: 19050
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
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