2

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.

4
  • Your example is missing the link. Jan 29, 2018 at 18:57
  • 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, 2018 at 19:39
  • @Luke Question reopened.
    – Aaron
    Jan 30, 2018 at 6:36
  • 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 at 14:53

2 Answers 2

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.

2
  • 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! Feb 2, 2018 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: stackoverflow.com/questions/7247298/size-of-list-in-memory Feb 2, 2018 at 16:53
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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