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'saffine_transform功能做到这一点？

4

1 回答 1

1

有几个选项。并不是所有的空间变换都在线性空间中，所以它们不能都使用仿射变换，所以不要总是依赖它。如果您有两个 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.affinity.affine_transform这与 Shapely函数中使用的基本算法相同。

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函数适用于图像或栅格数据，而不是矢量数据。

于 2013-04-15T23:32:39.457 回答