affine_transform xy координаты от gda94

Я пытаюсь понять, как преобразовать многоугольник, координаты которого находятся в Spatial Reference 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]]

Есть ли способ сделать это с помощью функции affine_transform scipy.ndimage?


python numpy scipy gdal
person James Mills    schedule 15.04.2013    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Есть несколько вариантов. Не все пространственные преобразования выполняются в линейном пространстве, поэтому не все они могут использовать аффинное преобразование, поэтому не всегда полагайтесь на него. Если у вас есть два EPSG SRID, вы можете выполнить общее пространственное преобразование с помощью модуля OSR GDAL. Некоторое время назад я написал пример, который можно адаптировать.

В противном случае аффинное преобразование имеет базовую математику:

                    / 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 предназначена для изображений или растровых данных, а не для векторных данных.

person Mike T    schedule 15.04.2013
comment
Да, модуль shapely.affinity отлично это реализует, и я могу использовать его для преобразования полигонов в EPSG 4283 (GDA94) в координаты xy для индексации в массив numpy (GeoTiff) (также в GDA94). На самом деле это не отвечает на вопрос, можно ли это сделать, преобразовав многоугольник в массив numpy и выполнив для него affine_transform напрямую (избегая циклов for и чистого python, который может быть медленнее для полигонов, содержащих › 500 000 точек). - person James Mills; 16.04.2013
comment
Я думаю, что мы с коллегой решили, что scipy.ndimage.affine_transform предназначен для работы с изображениями, а не с векторами/полигонами. Лучший способ достичь общего решения - выполнить affine_transform, например, как ваше решение выше; который реализует shapely.affinity.affine_transform. Я приму ваш ответ, если вы измените его, чтобы сказать, что это невозможно сделать с помощью scipy.ndimage и что вы должны выполнить аффинное преобразование вектора/полигона (аналогично вашему) или использовать shapely.affinity.affine_transform --James - person James Mills; 16.04.2013
comment
Спасибо, что помогли ответить на это :) - person James Mills; 16.04.2013