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中的光柵提取點

回答

相關問題