4

I am new to GDAL and am looking for pointers about the Correct™ way to accomplish a few tasks. I would like to do this in Python, not on the command line.

I have a GeoTIFF that shows elevation at 10m intervals over a given area. That is, every pixel is a 10m x 10m area. The GeoTIFF contains a single raster image, where elevation in meters is expressed as a Float32. The GeoTIFF seems to have tags defining the geoid/datum (which UTM zone etc).

I wish to open that image, and query what the elevation is at a certain point, such as in the following pseudocode:

elevation = myGeoTiffData.getValueAtPosition(Position('WGS84', 39.584588, -122.461320))

Further, I would like to have this elevation interpolated (linearly, spline-wise, however else) if the requested position does not fall on a 10m x 10m lattice point. I see there are some functions to do this (GDAL grid?), but I do not understand how to use them in this case.

I have tried to read the GDAL documentation, but I can't seem to figure this out. Is there any advice people here can offer?

3 Answers 3

7

A high-level rasterstats tool is built on GDAL to fetch data from raster and vector data sources.

import rasterstats
from shapely.geometry import Point

raster = "elevation.tif"
pt = Point(245000, 1000000)

# default is interpolate="bilinear"
elev = rasterstats.point_query([pt], raster)[0]

# or use nearest
elev = rasterstats.point_query([pt], raster, interpolate="nearest")[0]

Transformations to/from lat/long coordinates can be done using pyproj.

5

rioxarray has an interpolation method that will work really well for this. See the details here https://corteva.github.io/rioxarray/stable/examples/interpolate_na.html

import rioxarray

raster_file = './elevation.tif'

ds = rioxarray.open_rasterio(raster_file)
point_elevation = ds.interp(x= -122.461320, y= 39.584588, method='linear')
2
  • Is xarray aware of coordinate systems? Does it read the crs from the GeoTIFF tags?
    – user30184
    Commented Oct 28, 2022 at 6:29
  • 1
    @user30184 Yeah, it looks like xarray.open_rasterio and the newer rioxarray.open_rasterio use the rasterio library which does a bunch of nice GISsy things. Commented Oct 28, 2022 at 23:00
2

gdallocationinfo is a utility for getting a value from a raster band or bands given a point location. However it does not do interpolation. I would recommend resampling to a higher resolution with an appropriate resampling algorithm. gdalwarp is the standard GDAL tool for this, and includes several resampling algorithms but it is up to you to choose the appropriate one.

I'm unaware of any GDAL tool that can do an on-the-fly continuous point interpolation.

1
  • Thanks. I'm trying to do this all in Python, so I can do various computations on the data. I will update my question to reflect that Commented Oct 28, 2022 at 0:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.