python - 使用python和gdal从.geotiff中的像素点查找纬度/经度坐标
问题描述
我在 .geotiff imagesat 中有由像素坐标定义的点。
我正在尝试将这些像素坐标转换为经度和纬度。
起初,a 试图从我之前在这里提出的反问题中调整给定的代码
def pixel_to_world(geo_matrix, x, y):
ul_x = geo_matrix[0]
ul_y = geo_matrix[3]
x_dist = geo_matrix[1]
y_dist = geo_matrix[5]
_x = x * x_dist + ul_x
_y = y * y_dist + ul_y
return (_x, _y)
def build_transform_inverse(dataset, EPSG):
source = osr.SpatialReference(wkt=dataset.GetProjection())
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
return osr.CoordinateTransformation(source, target)
def find_spatial_coordinate_from_pixel(dataset, transform, x, y):
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(float(x), float(y))
point.Transform(transform)
return pixel_to_world(dataset.GetGeoTransform(), point.GetX(), point.GetY())
ds = gdal.Open(source_directory_path + filename)
_t = build_transform_inverse(ds, 4326)
coordinates = find_spatial_coordinate_from_pixel(ds, _t, point[0], point[1])
但我的最终坐标不正确,我得到类似(-1528281.9351183572, 1065990.778732022)而不是(-13.728790283203121 , 9.531686027524676)
ul_x = -1528281.943533814,
ul_y = 1065990.7964650677,
x_dist = 0.5970709765586435,
y_dist = -0.5972870511184641
我还尝试了光栅方法,例如:
import rasterio
with rasterio.open('path/to/file.tiff') as map_layer:
pixels2coords = map_layer.xy(2679,2157)
print(pixels2coords)
Witch 未能返回正确的长纬度坐标对(它返回 (-1526993.7629018887, 1064390.365811596) )。
我究竟做错了什么 ?
解决方案
您已经颠倒了在像素、投影坐标和经度之间转换的顺序。您正在按照pixels -> lat lngs -> projected
应有的顺序进行操作pixels -> projected -> lat lngs
。
以下应该工作
from osgeo import osr, ogr, gdal
def pixel_to_world(geo_matrix, x, y):
ul_x = geo_matrix[0]
ul_y = geo_matrix[3]
x_dist = geo_matrix[1]
y_dist = geo_matrix[5]
_x = x * x_dist + ul_x
_y = y * y_dist + ul_y
return _x, _y
def build_transform_inverse(dataset, EPSG):
source = osr.SpatialReference(wkt=dataset.GetProjection())
target = osr.SpatialReference()
target.ImportFromEPSG(EPSG)
return osr.CoordinateTransformation(source, target)
def find_spatial_coordinate_from_pixel(dataset, transform, x, y):
world_x, world_y = pixel_to_world(dataset.GetGeoTransform(), x, y)
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(world_x, world_y)
point.Transform(transform)
return point.GetX(), point.GetY()
ds = gdal.Open(source_directory_path + filename)
_t = build_transform_inverse(ds, 4326)
coordinates = find_spatial_coordinate_from_pixel(ds, _t, point[0], point[1])
print(coordinates)
在find_spatial_coordinate_from_pixel
我已经翻转了顺序,pixel_to_world
现在在坐标变换之前首先调用它。
推荐阅读
- tensorflow - DLL 加载失败:找不到指定的过程。并且没有显示任何丢失的 DLL 文件名
- python-3.x - 如何使用“while”语句不断地向用户提出相同的问题?
- python - How to parse the heatmap output for the pose estimation tflite model?
- javascript - 您如何使用 CSS 将三个子 div 居中在父 div 中?
- node.js - NodeJS xml-stream,解析未知节点
- java - 在 java 的继承中以不同的方式创建对象时会发生什么?
- encryption - Openssl 1.0.2p 解密失败?
- android - Android:我已经删除了我的数据,但房间里仍然存在一些数据
- c++ - 未知的 CMake 命令“create_single_source_cgal_program”
- sql - 在 SQL Server 2012 中每次都询问未提交的事务