Extract raster values from point using GDAL

Question

I want to extract raster value from point shapefile geometry. Raster size is 45924 x 61671 and shapefile has 11.000.000 points. I'm using extract function from Raster package in R but is not very fast.

I'd like avoid non-terminal methods like ArcGIS or QGIS and prefer some terminal options like gdal_xx.py mode at prompt.

As a result I need an object with 11.000.000 length in any format (CSV, txt, shp, etc).

I need repeat this procedure over several similar rasters.

Here's a nice example but I don't know how to get or manipulate the resulting object and I'm not sure about their efficiency.

There appears to be two questions here--one a duplicate and another a follow-up to a previous question: 1) how to manipulate the resulting object from a GDAL sample operation? and 2) How to sample rasters with OGR point? Please edit your question to clarify. — Aaron, Jan 29 '18 at 19:39
I'd dump XY point coordinates and use gdallocation to get the raster values based on the coordinates. If you put all of that in a for loop as a bash/shell script (or directly in the terminal) you can add as many files as you want. — Gery, Mar 23 '22 at 14:53

score 5 · Answer 1 · edited Jun 09 '18 at 06:55

5

Working from your provided example (which is the standard way of doing this), you can collect all the info in a python list.

Let's say your point layer has a unique ID field (if it doesn't, create one, as it really should). For this example, let's call it "id_points". You complement the code in your link with:

li_values = list()
for feat in lyr:
    geom = feat.GetGeometryRef()
    feat_id = feat.GetField('id_points')
    mx, my = geom.GetX(), geom.GetY()

    px = int((mx - gt[0]) / gt[1])
    py = int((my - gt[3]) / gt[5])

    intval = rb.ReadAsArray(px, py, 1, 1)
    li_values.append([feat_id, intval[0]])

Original code attribution

This gives you a list of feature IDs and their associated raster values. You can then save it in a CSV (for example):

import csv

with open(r'csv/file/path.csv', 'wb') as csvfile:
    wr = csv.writer(csvfile)
    wr.writerows(li_values)

This will give you an output in the form of a table, which you can then open anywhere you want.

edited Jun 09 '18 at 06:55

user2856

65,736
6
115
196

answered Jan 29 '18 at 19:27

Roberto Ribeiro

3,161
13
28

Thks, Roberto. I have a little concern about the line li_values.append([feat_id, intval[0]]). Its possible to create an empty object with length X before the loop starts and fill it in each iteration? I dont know if in each iteration the PC have to create an object bigger and bigger under this scheme. In R (sorry but Im not so good in Python), for example reesults faster "li_values[positionX] <- intval[0]" than "li_values <- c(li_values, intval[0)]" or "li_values <- append(li_values, intval[0)]" (c, append, rbind, cbind, etc). Thanks! – gonzalez.ivan90 Feb 02 '18 at 16:33
@gonzalez.ivan90 You can using numpy arrays instead of lists, but they perform the same. Python list appending has very little overhead, you can use it without fear. More info here: https://stackoverflow.com/questions/7247298/size-of-list-in-memory – Roberto Ribeiro Feb 02 '18 at 16:53

score 0 · Answer 2 · edited Mar 23 '22 at 14:33

Using multiple solutions from other users I've come up with this:

from math import floor
from osgeo import gdal,ogr
import csv
src_filename = '/path of your tif file.tif'
shp_filename = '/path of your shape file.shp'
src_ds=gdal.Open(src_filename) 
gt_forward=src_ds.GetGeoTransform()
gt_reverse=gdal.InvGeoTransform(gt_forward)
rb=src_ds.GetRasterBand(1)
ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
feat = lyr.GetFeature(1)
li_values = list()
for feat in lyr:
  geom = feat.GetGeometryRef()
  feat_id = feat.GetField('id_points')
  mx, my = geom.GetX(), geom.GetY()
px = int((mx - gt_forward[0]) / gt_forward[1])
  py = int((my - gt_forward[3]) / gt_forward[5])
intval = rb.ReadAsArray(px, py, 1, 1)
  li_values.append([feat_id, intval[0]])
with open(r'/path.csv', 'w') as csvfile:
    wr = csv.writer(csvfile)
    wr.writerows(li_values)

this doesn't work properly – Flash Thunder Jun 29 '22 at 09:18 — Flash Thunder, Jun 29 '22 at 09:18

Extract raster values from point using GDAL

2 Answers2

Linked