从GDAL中的栅格中提取点

Posted 2023-02-23

【中文标题】从GDAL中的栅格中提取点

【英文标题】：Extract Point From Raster in GDAL

【发布时间】：2012-11-18 10:35:52

【问题描述】：

我有一个光栅文件和一个 WGS84 纬度/经度点。

我想知道栅格中的哪个值与该点对应。

我的感觉是我应该在光栅对象或其一个波段上使用GetSpatialRef()，然后将ogr.osr.CoordinateTransformation() 应用于该点以将其映射到光栅空间。

然后我希望我可以简单地询问光栅乐队当时的情况。

但是，光栅对象似乎没有GetSpatialRef() 或访问地理定位点的方法，所以我有点不知所措。

有什么想法吗？

【参考方案1】：

假设我有一个 geotiff 文件 test.tif。然后下面的代码应该在像素附近的某个地方查找值。我对查找单元格的部分没有那么自信，并且会修复存在错误。这个页面应该有帮助，"GDAL Data Model"

另外，如果你还没有，你可以去gis.stackexchange.com找专家。

import gdal, osr

class looker(object):
    """let you look up pixel value"""

    def __init__(self, tifname='test.tif'):
       """Give name of tif file (or other raster data?)"""

        # open the raster and its spatial reference
        self.ds = gdal.Open(tifname)
        srRaster = osr.SpatialReference(self.ds.GetProjection())

        # get the WGS84 spatial reference
        srPoint = osr.SpatialReference()
        srPoint.ImportFromEPSG(4326) # WGS84

        # coordinate transformation
        self.ct = osr.CoordinateTransformation(srPoint, srRaster)

        # geotranformation and its inverse
        gt = self.ds.GetGeoTransform()
        dev = (gt[1]*gt[5] - gt[2]*gt[4])
        gtinv = ( gt[0] , gt[5]/dev, -gt[2]/dev, 
                gt[3], -gt[4]/dev, gt[1]/dev)
        self.gt = gt
        self.gtinv = gtinv

        # band as array
        b = self.ds.GetRasterBand(1)
        self.arr = b.ReadAsArray()

    def lookup(self, lon, lat):
        """look up value at lon, lat"""

        # get coordinate of the raster
        xgeo,ygeo,zgeo = self.ct.TransformPoint(lon, lat, 0)

        # convert it to pixel/line on band
        u = xgeo - self.gtinv[0]
        v = ygeo - self.gtinv[3]
        # FIXME this int() is probably bad idea, there should be 
        # half cell size thing needed
        xpix =  int(self.gtinv[1] * u + self.gtinv[2] * v)
        ylin = int(self.gtinv[4] * u + self.gtinv[5] * v)

        # look the value up
        return self.arr[ylin,xpix]

# test
l = looker('test.tif')
lon,lat = -100,30
print l.lookup(lon,lat)

lat,lon =28.816944, -96.993333
print l.lookup(lon,lat)

【讨论】：

我不遵循用于在构造函数中计算 dev 或 gtinv 的逻辑。你可以解释吗？此外，很高兴知道如何在您的 ImportFromEPSG(4326) 语句中处理不同的预测。谢谢。

我相信我有地理变换矩阵的逆矩阵，然后当我得到 lon/lat 时，我反算了像素/线

【参考方案2】：

是的，API 不一致。栅格（数据源）有一个 GetProjection() 方法（返回 WKT）。

这是一个做你想做的函数（取自here）：

def extract_point_from_raster(point, data_source, band_number=1):
    """Return floating-point value that corresponds to given point."""

    # Convert point co-ordinates so that they are in same projection as raster
    point_sr = point.GetSpatialReference()
    raster_sr = osr.SpatialReference()
    raster_sr.ImportFromWkt(data_source.GetProjection())
    transform = osr.CoordinateTransformation(point_sr, raster_sr)
    point.Transform(transform)

    # Convert geographic co-ordinates to pixel co-ordinates
    x, y = point.GetX(), point.GetY()
    forward_transform = Affine.from_gdal(*data_source.GetGeoTransform())
    reverse_transform = ~forward_transform
    px, py = reverse_transform * (x, y)
    px, py = int(px + 0.5), int(py + 0.5)

    # Extract pixel value
    band = data_source.GetRasterBand(band_number)
    structval = band.ReadRaster(px, py, 1, 1, buf_type=gdal.GDT_Float32)
    result = struct.unpack('f', structval)[0]
    if result == band.GetNoDataValue():
        result = float('nan')
    return result

其文档如下（取自here）：

spatial.extract_point_from_raster(point, data_source, band_number=1)

data_source 是 GDAL 栅格，point 是 OGR 点对象。这 函数返回指定波段的像素值 离点最近的data_source。

point 和 data_source 不必在同一个参考系中，但 它们都必须定义适当的空间参考。

如果点不在光栅中，则引发 RuntimeError。

【参考方案3】：

project = self.ds.GetProjection()
srPoint = osr.SpatialReference(wkt=project)

完成...这样，矢量文件采用了输入光栅文件的投影