GDAL中的光栅提取点

我有一个光栅文件和一个WGS84经纬度点。GDAL中的光栅提取点

我想知道栅格中的值与点的对应关系。

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

我的希望是，那么我可以简单地问一下栅格乐队的那一点。

但是，栅格对象似乎没有GetSpatialRef()或访问位于地理位置的点的方式，所以我对如何执行此操作有些遗憾。

有什么想法？

来源

2012-11-18 Richard

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

来源

2012-11-19 00:21:38 yosukesabai

我不遵循用于在构造函数中计算'dev'或'gtinv'的逻辑。你可以解释吗？另外，在'ImportFromEPSG（4326）'语句中知道如何处理不同的投影会很好。谢谢。 – theJollySin

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

-1

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

做......就这样，矢量文件已经通过从输入光栅文件

来源

2013-11-26 14:45:42 Emiliano

投影是的，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栅格和点是OGR点对象。函数返回距离点最近的指定波段的数据源的像素值。

point和data_source不需要在同一参考系统中，但它们必须都具有已定义的适当空间参考。

如果点不落在栅格中，则会引发RuntimeError。

来源

2015-01-30 14:30:21

GDAL中的光栅提取点

回答

相关问题