python - GIS/GEOTiff/GDAL/Python 如何从像素获取坐标

Question

我正在研究从 GEOTiff 文件中检测对象并返回对象坐标的项目，这些输出将用于无人机飞向这些坐标

我使用带有 YOLO v2(图像检测器框架)和 OpenCV 的 tensorflow 来检测我在 GEOTiff 中需要的对象

import cv2
from darkflow.net.build import TFNet
import math
import gdal

# initial stage for YOLO v2 
options = {
    'model': 'cfg/yolo.cfg',
    'load': 'bin/yolov2.weights',
    'threshold': 0.4,
}
tfnet = TFNet(options)

# OpenCV read Image
img = cv2.imread('final.tif', cv2.IMREAD_COLOR)
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)

#Predict the image
result = tfnet.return_predict(img)

#Calculate the center and radius of each objects
i = 0
while i < len(result):
    tl = (result[i]['topleft']['x'], result[i]['topleft']['y'])
    br = (result[i]['bottomright']['x'], result[i]['bottomright']['y'])
    point = (int((result[i]['topleft']['x']+result[i]['bottomright']['x'])/2), int((result[i]['topleft']['y']+result[i]['bottomright']['y'])/2))
    radius = int(math.hypot(result[i]['topleft']['x'] - point[0], result[i]['topleft']['y'] - point[1]))
    label = result[i]['label']
    result[i]['pointx'] = point[0]
    result[i]['pointy'] = point[1]
    result[i]['radius'] = radius    
    i += 1

print(result)

所以结果出来就像一组JSON

[{'label': 'person', 'confidence': 0.6090355, 'topleft': {'x': 3711, 'y': 1310}, 'bottomright': {'x': 3981, 'y': 1719}, 'pointx': 3846, 'pointy': 1514, 'radius': 244}]

如您所见，对象的位置以像素 (x,y) 为单位返回我想用这些 x,y 转换为 lat,lng 中的坐标所以我尝试使用 GDAL(用于读取图像中包含的 GEO 信息的库)

所以这是在终端中使用 gdalinfo 的图像的 GEO 信息

Driver: GTiff/GeoTIFF
Files: final.tif
Size is 8916, 6888
Coordinate System is:
PROJCS["WGS 84 / UTM zone 47N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32647"]]
Origin = (667759.259870000067167,1546341.352920000208542)
Pixel Size = (0.032920000000000,-0.032920000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  667759.260, 1546341.353) (100d33'11.42"E, 13d58'57.03"N)
Lower Left  (  667759.260, 1546114.600) (100d33'11.37"E, 13d58'49.65"N)
Upper Right (  668052.775, 1546341.353) (100d33'21.20"E, 13d58'56.97"N)
Lower Right (  668052.775, 1546114.600) (100d33'21.15"E, 13d58'49.59"N)
Center      (  667906.017, 1546227.976) (100d33'16.29"E, 13d58'53.31"N)
Band 1 Block=8916x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
Band 2 Block=8916x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
Band 3 Block=8916x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
Band 4 Block=8916x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

有人吗？

Answer 1

您需要使用与光栅文件关联的 GeoTransform 矩阵将像素坐标转换为地理空间。使用 GDAL，您可以执行以下操作:

# open the dataset and get the geo transform matrix
ds = gdal.Open('final.tif') 
xoffset, px_w, rot1, yoffset, px_h, rot2 = ds.GetGeoTransform()

# supposing x and y are your pixel coordinate this 
# is how to get the coordinate in space.
posX = px_w * x + rot1 * y + xoffset
posY = rot2 * x + px_h * y + yoffset

# shift to the center of the pixel
posX += px_w / 2.0
posY += px_h / 2.0

当然，您获得的位置将相对于用于栅格数据集的同一坐标引用系。因此，如果您需要将其转换为纬度/经度，则必须做进一步的阐述:

# get CRS from dataset 
crs = osr.SpatialReference()
crs.ImportFromWkt(ds.GetProjectionRef())
# create lat/long crs with WGS84 datum
crsGeo = osr.SpatialReference()
crsGeo.ImportFromEPSG(4326) # 4326 is the EPSG id of lat/long crs 
t = osr.CoordinateTransformation(crs, crsGeo)
(lat, long, z) = t.TransformPoint(posX, posY)

抱歉，我对 Python 不是很流利，所以您可能需要调整此代码。查看 GeoTransform 的文档 here for the C++ API了解有关矩阵元素的更多信息。