从GDAL中的光栅提取点

Question

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

我想知道光栅中的值与点对应。

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

我的希望是我可以简单地问拉斯特乐队在那一点上是什么。

然而，光栅对象似乎没有GetSpatialRef()或访问地理位置点的方法，所以我对如何做到这一点有些茫然。

有什么想法吗？

Answer 1

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

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

Answer 2

假设我有一个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)

Answer 3

是的，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 is a GDAL raster, and point is an OGR point object. The function returns the value of the pixel of the specified band of data_source that is nearest to point.

point and data_source need not be in the same reference system, but they must both have an appropriate spatial reference defined.

If the point does not fall in the raster, RuntimeError is raised.