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.

  • Your example is missing the link. – Roberto Ribeiro Jan 29 '18 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 '18 at 19:39
  • @Luke Question reopened. – Aaron Jan 30 '18 at 6:36

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)

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

  • 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 2 '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: stackoverflow.com/questions/7247298/size-of-list-in-memory – Roberto Ribeiro Feb 2 '18 at 16:53

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.