Sample raster cell value at coordinates?

Question

I have many rasters in ESRI ASCII format (http://resources.arcgis.com/en/help/main/10.1/index.html#/Esri_ASCII_raster_format/009t0000000z000000/).

I need to extract cell values at given locations / coordinates. Can anyone suggest a python package to achieve this? I suspect there may be something in gdal tools but I have been unable to find anything so far.

I am looking for similar functionality to GMT grdtrack, with which you can pass a table of coordinates and retrieve the cell value. http://gmt.soest.hawaii.edu/doc/5.1.0/grdtrack.html

However I was hoping / wondering if this is possible in python as my previous and later stages of my analysis are all in python.

Questions asking us to recommend or find a book, tool, software library, tutorial or other off-site resource are off-topic for Stack Overflow as they tend to attract opinionated answers and spam. Instead, describe the problem and what has been done so far to solve it. — , Nov 08 '15 at 16:52

Constantinius · Accepted Answer · 2020-03-21T08:29:03.560

1

Of course, you can do this with GDAL. Provided that you have the coordinate where you want to sample your raster in the same projection you can do something like this:

from osgeo import gdal

def world2Pixel(gt, x, y):
  ulX = gt[0]
  ulY = gt[3]
  xDist = gt[1]
  yDist = gt[5]
  rtnX = gt[2]
  rtnY = gt[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / yDist)
  return (pixel, line)

dataset = gdal.Open(filename)
gt = dataset.GetGeoTransform()

pixel, line = world2Pixel(gt, x, y)

band = dataset.GetRasterBand(1)
value = band.ReadAsArray(pixel, line, 1, 1)[0, 0]

edited Mar 21 '20 at 08:29

answered Nov 08 '15 at 17:05

Constantinius

34,183
8
77
85

Thanks that works great! Apologies to all for posting an ambiguous question: my original solution used GMT grdtrack. "grdtrack $table_of_XY -G$grid -Ql > $values_out" – BenP Nov 09 '15 at 18:13
There seems to be an error in the answer, please verify if this: `line = int((ulY - y) / xDist)` should be this: `line = int((ulY - y) / yDist)` – Argantonio65 Mar 19 '20 at 12:08
It should actually be: `line = int((ulY - y) / - yDist)`, since gt[5] is returned as a negative value (as it reads from here https://gdal.org/tutorials/raster_api_tut.html). Sorry for the double update. Actually, I saw that here: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html they use the same world2Pixel function in the incorrect form. – Argantonio65 Mar 22 '20 at 17:19

Sample raster cell value at coordinates?

1 Answers1