在坐标处采样栅格像元值？

Question

我有很多 ESRI ASCII 格式的栅格 (http://resources.arcgis.com/en/help/main/10.1/index.html#/Esri_ASCII_raster_format/009t0000000z000000/)。

我需要提取给定位置/坐标处的单元格值。任何人都可以建议 python 包来实现这一目标吗？我怀疑 gdal 工具中可能有一些东西，但到目前为止我找不到任何东西。

我正在寻找与 GMT grdtrack 类似的功能，您可以使用它传递 table 坐标并检索单元格值。 http://gmt.soest.hawaii.edu/doc/5.1.0/grdtrack.html

但是我希望/想知道这在 python 中是否可行，因为我之前和后期的分析阶段都在 python.

中

Answer 1

当然，您可以使用 GDAL 来完成此操作。如果您拥有要在同一投影中对栅格进行采样的坐标，您可以执行以下操作：

from osgeo import gdal

def world2Pixel(gt, x, y):
  ulX = gt[0]
  ulY = gt[3]
  xDist = gt[1]
  yDist = gt[5]
  rtnX = gt[2]
  rtnY = gt[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / yDist)
  return (pixel, line)

dataset = gdal.Open(filename)
gt = dataset.GetGeoTransform()

pixel, line = world2Pixel(gt, x, y)

band = dataset.GetRasterBand(1)
value = band.ReadAsArray(pixel, line, 1, 1)[0, 0]

在坐标处采样栅格像元值？

Sample raster cell value at coordinates?

python

gdal

gmt