python - 从 GDAL 中的栅格中提取点

Question

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

我想知道栅格中的什么值与该点相对应。

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

我的希望是我可以简单地询问光栅的波段此时是什么。

但是，栅格对象似乎没有 GetSpatialRef() 或访问地理定位点的方法，因此我有点不知道如何执行此操作。

有什么想法吗？

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