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

Question

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

我使用 tensorflow 与 YOLO v2（图像检测器框架）和 OpenCV 来检测我在 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 的文档以了解有关矩阵元素的更多信息。

Answer 2

如果没有 Gabriella 发布的出色而清晰的 Py​​thon 代码，我不知道我是否会想出如何在 C 中做到这一点。gdal 的文档和示例非常稀少。

这是 Gabriella 代码的 C 版本：

const char fn[] = "/path/to/geo/file.tif";

GDALDatasetH  hDataset;
GDALAllRegister(); // Register all GDAL formats
hDataset = GDALOpen( fn, GA_ReadOnly ); // Open our geo file (GeoTIFF or other supported format)
if (hDataset == NULL)
{
    printf("Failed to open dataset\n");
    return;
}

// These are the input points to be transformed, in pixel coordinates of the source raster file
double x = 20;
double y = 20;

double        adfGeoTransform[6];
GDALGetGeoTransform( hDataset, adfGeoTransform );

// Put the returned transform values into named vars for readability
double xoffset = adfGeoTransform[0];
double px_w = adfGeoTransform[1];
double rot1 = adfGeoTransform[2];
double yoffset = adfGeoTransform[3];
double rot2 = adfGeoTransform[4];
double px_h = adfGeoTransform[5];

// Apply transform to x,y. Put into posX,posY
double posX = px_w * x + rot1 * y + xoffset;
double posY = rot2 * x + px_h * y + yoffset;

// Transform to center of pixel
posX += px_w / 2.0;
posY += px_h / 2.0;

OGRErr err = 0;

// sr0 is the "from" spatial reference, pulled out of our file
OGRSpatialReferenceH sr0 = OSRNewSpatialReference(GDALGetProjectionRef(hDataset));

// sr1 is the "to" spatial reference, initialized as EPSG 4326 (lat/lon)
OGRSpatialReferenceH sr1 = OSRNewSpatialReference(NULL);
err = OSRImportFromEPSG(sr1, 4326);

double xtrans = posX;
double ytrans = posY;
double ztrans = 0;
int pabSuccess = 0;

// Make our transformation object
OGRCoordinateTransformationH trans = OCTNewCoordinateTransformation(sr0, sr1);

// Transform our point posX,posY, put it into xTrans,yTrans
OCTTransformEx(trans, 1, &xtrans, &ytrans, &ztrans, &pabSuccess);
    
GDALClose(hDataset);

printf("map coordinates (%f, %f)\n", xtrans, ytrans);