[GIS] Extract raster values from point using GDAL

extractgdalraster

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.

Best Answer

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.

Related Question