给定一个 geotiff 文件，如何找到最接近给定 latitude/longitude 的单个像素？

Question

我有一个要在 Python 中使用 gdal 打开的 geotiff 文件，我需要找到最接近指定 latitude/longitude 的单个像素。我之前使用的是类似数据的不相关文件类型，所以我对 gdal 和 geotiff 都是全新的。

如何做到这一点？到目前为止我所拥有的是

import gdal
ds = gdal.Open('foo.tiff')
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
gp = ds.GetProjection()
data = np.array(ds.ReadAsArray())
print(gt)
print(gp)

生成（我的文件）

(-3272421.457337171, 2539.703, 0.0, 3790842.1060354356, 0.0, -2539.703)

和

PROJCS["unnamed",GEOGCS["Coordinate System imported from GRIB file",DATUM["unnamed",SPHEROID["Sphere",6371200,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",25],PARAMETER["central_meridian",265],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",25],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

理想情况下，应该有一个简单的函数调用，它还会 return 指示指定位置是否落在栅格范围之外。

我的后备方案是从另一个来源获取包含每个像素的纬度和经度的网格，然后对所需位置进行蛮力搜索，但我希望有更优雅的方法。

注意：我认为我正在尝试做的相当于命令行

gdallocationinfo -wgs84 foo.tif <longitude> <latitude>

其中 return 的结果类似于

Report:
  Location: (1475P,1181L)
  Band 1:
    Value: 66

这向我表明，该功能可能已经在 gdal 模块中，如果我能找到正确的调用方法的话。

Answer 1

你基本上需要两个步骤：

将 lat/lon 点转换为光栅投影
将 mapx/mapy（在光栅投影中）转换为像素坐标

鉴于您已经在上面发布的代码，可以通过以下方式定义两个投影系统：

from osgeo import gdal, osr

point_srs = osr.SpatialReference()
point_srs.ImportFromEPSG(4326) # hardcode for lon/lat

# GDAL>=3: make sure it's x/y
# see https://trac.osgeo.org/gdal/wiki/rfc73_proj6_wkt2_srsbarn
point_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)     

file_srs = osr.SpatialReference()
file_srs.ImportFromWkt(gp)

创建坐标转换，并使用它将点从 lon/lat 转换为 mapx/mapy 坐标（无论投影是什么）：

ct = osr.CoordinateTransformation(point_srs, file_srs)

point_x = -114.06138 # lon
point_y = 51.03163  # lat
mapx, mapy, z = ct.TransformPoint(point_x, point_y)

要从地图坐标转换为像素坐标，需要先反转地理变换。然后可以用来检索像素坐标，如：

gt_inv = gdal.InvGeoTransform(gt)
pixel_x, pixel_y = gdal.ApplyGeoTransform(gt_inv, mapx, mapy)

四舍五入这些像素坐标应该允许您使用它们来索引数据数组。如果您查询的点在栅格之外，您可能需要裁剪它们。

# round to pixel
pixel_x = round(pixel_x)
pixel_y = round(pixel_y)

# clip to file extent
pixel_x = max(min(pixel_x, width-1), 0)
pixel_y = max(min(pixel_y, height-1), 0)

pixel_data = data[pixel_y, pixel_x]

给定一个 geotiff 文件，如何找到最接近给定 latitude/longitude 的单个像素？

Given a geotiff file, how does one find the single pixel closest to a given latitude/longitude?

python

gdal

geotiff