从GDAL中的栅格中提取点 [英] Extract Point From Raster in GDAL

问题描述

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

I have a raster file and a WGS84 lat/lon point.

我想知道栅格中的哪个值与该点相对应.

I would like to know what value in the raster corresponds with the point.

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

My feeling is that I should use GetSpatialRef() on the raster object or one of its bands and then apply a ogr.osr.CoordinateTransformation() to the point to map it to the raster's space.

那么我的希望是，我可以简单地问一下栅格乐队当时的情况.

My hope would then be that I could simply ask the rasters' bands what is at that point.

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

However, the raster object doesn't seem to have a GetSpatialRef() or a way to access a geo-located point, so I'm somewhat at a loss for how to do this.

有什么想法吗?

推荐答案

说我有一个geotiff文件test.tif.然后跟随代码应在像素附近的某个位置查找值.我对查找单元格的部分不太自信，并且会修复存在的错误.此页面应该有帮助，"GDAL数据模型"

Say i have a geotiff file test.tif. Then followin code should look up value somewhere near the pixel. I am not that confident for the part looking up cell, and will fix there is error. This page should help, "GDAL Data Model"

如果您还没有的话，也可以访问 gis.stackexchange.com 来找到专家.

Also, you may go to gis.stackexchange.com to find experts, if you haven't.

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)

这篇关于从GDAL中的栅格中提取点的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！