GDA94坐标系下的仿射变换xy坐标

Question

GDA94坐标系下的仿射变换xy坐标

3

我正在尝试弄清楚如何将坐标系为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函数有没有办法做到这一点？

- James Mills

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mike T · Accepted Answer

有几个选项。并非所有的空间转换都在线性空间中，因此它们不能全部使用仿射变换，因此不要总是依赖它。如果您有两个EPSG SRID，则可以使用GDAL的OSR模块进行通用空间转换。我以前写过一个示例，可以进行调整。

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

                    / 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函数相同的基本算法。

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函数旨在处理图像或栅格数据，而不是矢量数据。