我有一个光栅文件和一个WGS84纬度/经度点。
我想知道栅格中的哪个值与该点相对应。
我的感觉是我应该在栅格对象或其波段之一上使用GetSpatialRef()
,然后在该点上应用ogr.osr.CoordinateTransformation()
以将其映射到栅格空间。
我的希望是,我可以简单地问一下栅格乐队当时的状况。
但是,栅格对象似乎没有GetSpatialRef()
或访问地理定位点的方法,因此我对如何执行此操作有些困惑。
有什么想法吗?
最佳答案
说我有一个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)
关于python - 从GDAL中的栅格提取点,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/13439357/