来自 gda94 的 affine_transform xy 坐标

Question

我试图弄清楚如何将坐标位于空间参考 GDA94 (EPSG 4283) 中的多边形转换为 xy 坐标（反仿射变换矩阵）。

以下代码有效：

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

这是打印矩阵的输出：

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

有没有办法使用scipy.ndimage 的affine_transform 函数来做到这一点？

Answer 1

有几个选项。并不是所有的空间变换都在线性空间中，所以它们不能都使用仿射变换，所以不要总是依赖它。如果您有两个 EPSG SRID，您可以使用 GDAL 的 OSR 模块进行通用空间变换。 I wrote an example a while back，可适配。

---

否则，仿射变换具有基本数学：

                    / 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

可以在 Python 中通过点列表实现。

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

这与 Shapely 的 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)

最后，值得一提的是，scipy.ndimage.interpolation.affine_transform function 用于图像或栅格数据，而不是矢量数据。