GIS/GEOTiff/GDAL/Python如何从像素获取坐标 [英] GIS / GEOTiff / GDAL / Python How to get coordinates from pixel

问题描述

我正在研究从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

有人吗?

解决方案

您需要使用与栅格文件关联的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的文档此处是C ++ API ，以了解有关矩阵元素的更多信息. /p>

I'm working on the project to detect the object from GEOTiff files and return coordinates of the objects and those output will use for drone to fly to those coordinate

I use tensorflow with YOLO v2(image detector framework) and OpenCV to detect the objects that I need in 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)

So the results come out like set of JSON

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

as you can see the location of the object is return in pixel (x,y) and I want to use these x,y to convert to coordinate in lat,lng so I try to use GDAL (the library use for read the GEO infomation that contain inside the image)

so here's the GEO infomation of the image by using gdalinfo in terminal

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

Any one?

解决方案

You need to transform pixel coordinates to geographic space using the GeoTransform matrix that is associated to your raster files. Using GDAL you could do something like the following:

# 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

Of course the position you get will be relative to the same coordinate reference system that is used for your raster dataset. So if you need to transform it to lat/long you will have to do further elaborations:

# 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)

Sorry I'm not really fluent in python, so probably you will have to adapt this code. Checkout the documentation of GeoTransform here for the C++ API to learn more about the matrix elements.

这篇关于GIS/GEOTiff/GDAL/Python如何从像素获取坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！