Extract point from raster in GDAL

Question

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

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

It seems to me that I should use GetSpatialRef() for the raster object or one of its groups, and then apply ogr.osr.CoordinateTransformation() to the point to map it to the raster space.

My hope would then be that I could just ask the raster groups that at this point.

However, the raster object does not have GetSpatialRef() or a way to access a point in geolocation, so I don’t understand a bit how to do this.

Any thoughts?

Answer 1

Let's say I have a geotiff file test.tif. Then the followin code should look for a value somewhere near the pixel. I'm not sure if the part looking at the cell will fix the error. This page should help, "GDAL Data Model"

Alternatively, you can 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) 

Answer 2

Yes, the API is incompatible. The raster (data source) has instead the GetProjection() method (which returns WKT).

Here is the function that does what you want (picture 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 

Its documentation is as follows (figure here ):

 spatial.extract_point_from_raster(point, data_source, band_number=1) 

data_source is a GDAL raster, and a dot is an OGR object. the function returns the pixel value of the specified range data_source, closest to the point.

the point and the data source must not be in the same reference frame, but they must have the appropriate spatial reference.

If the point does not fall into the raster, a RuntimeError raises.

Answer 3

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

done ... with this, the vector file accepted the projection from the input raster file