GIS / GEOTiff / GDAL / Python 如何从像素获取坐标
GIS / GEOTiff / GDAL / Python How to get coordinates from pixel
我正在研究从 GEOTiff 文件和 return 对象坐标中检测对象的项目，这些输出将用于无人机飞向这些坐标
我使用带有 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}]
如您所见，对象的位置是 return 像素 (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
有人吗？
您需要使用与光栅文件关联的 GeoTransform 矩阵将像素坐标转换为地理坐标 space。使用 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
当然，您获得的位置将相对于用于栅格数据集的同一坐标参考系统。因此，如果您需要将其转换为 lat/long，您将需要做进一步的阐述：
# 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 以了解有关矩阵元素的更多信息。
如果没有 Gabriella 发布的出色而清晰的 Python 代码，我不知道我是否会想出如何在 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);
我正在研究从 GEOTiff 文件和 return 对象坐标中检测对象的项目，这些输出将用于无人机飞向这些坐标
我使用带有 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}]
如您所见，对象的位置是 return 像素 (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
有人吗？
您需要使用与光栅文件关联的 GeoTransform 矩阵将像素坐标转换为地理坐标 space。使用 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
当然，您获得的位置将相对于用于栅格数据集的同一坐标参考系统。因此，如果您需要将其转换为 lat/long，您将需要做进一步的阐述：
# 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 以了解有关矩阵元素的更多信息。
如果没有 Gabriella 发布的出色而清晰的 Python 代码，我不知道我是否会想出如何在 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);