Extracting raster values at points using Open Source GIS?

Question

How can I extract values from a raster by points?

I prefer not in Arcgis.

I prefer in Qgis or Mapwindow or other open source gis.

Answer 1

QGIS "Point Sampling Tool" should be the plugin you're looking for.

Here's a detailed description of how to use it: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Update based on Paolo's comment:

the plugin is not the only solution, and not always the easiest solution anymore. An alternative solution is the Saga function ‘Add raster values to point’ in the processing toolbox. See for details http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

Answer 2

I was having problems with the QGIS and SAGA GUI tools mentioned in this thread (Raster values to points was failing for some reason and throwing unhelpful errors and the GRASS v.sample created a whole new layer which was not helpful). After failing with the GUI tools for a while, I tried doing this in the Field Calculator. It worked quite well and I was able to control the process a little better than the GUIs allow, and make some other calculations along the way.

Say you have a layer named pts and another named rast, both in the same coordinate system. You'd like to sample rast at each X,Y pair represented in pts.

If you haven't used the Field Calculator before, it's pretty simple. You will enter your calculation in the "Expression" box, and Q gives you a number of variables and operations in the middle column, with help text on in the right column. I'll break this process into four steps:

1. Open the attribute table of the pts layer you'd like to sample with.

2. Once you are in the Field Calculator dialog, choose whether you'd like to Create a new field or Modify an existing field in your pts layer.

3. Next, build an expression to fill the new or existing pts attribute column. You might start by modifying the expression code that worked for me:

raster_value('rast', 1, make_point($x, $y))

4. You must supply raster_value() with a raster layer name 'rast', a band number 1, and the point geometry at make_point(). $x and $y are geometry variables dependent on the location of the point in each row of the attribute table.

This method also allows arithmetic operations like subtracting the value of another raster layer called other_rast from rast, which saved me a bunch of time over the GUI tools. Example below:

raster_value('rast', 1, make_point($x, $y)) - raster_value('other_rast', 1, make_point($x, $y))

Note again that the three layers pts, rast, and other_rast must be in the same coordinate system for this method to work.

Answer 3

In PostGIS 2.0 you can do:

SELECT ST_Value(rast, geom) val
FROM yourrastertabe, yourpointtable
WHERE ST_Intersects(rast, geom)

Make sure your raster is tiled very small when you load it (-t 10x10 with the loader).

Answer 4

Try using QGIS 3.2.2 and SAGA (installed by default in QGIS): "Raster Values to Points" function will do everything for you: It takes a image file and converts it into a Point-vector shape taking the information from raster image.

Answer 5

Hawthorne Beyer's GME tools do this nicely via command line, and allow easy batching with 'for' loops.

isectpntrst(in="path/to/shapefile", raster="path/to/raster", field="fieldname")

GME isectpntrst command reference

Answer 6

In GRASS GIS, you can either query the map in the GUI or use http://grass.osgeo.org/gdp/html_grass64/r.what.html

Answer 7

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html has step by step guide to use R Raster package extract raster values from points.

Answer 8

You can use this: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

Its in the SAGA Toolbox of Qgis! It does everything in one step :)

Answer 9

Here's a function I wrote using python and gdal. The function takes a list of rasters and a pandas dataframe containing the point coordinates and returns a pandas dataframe with the point coordinates, the centroids for the respective raster cells and the respective cell values. The function is part of the under development package chorospy package (found here).

import pandas
import numpy
from osgeo import gdal

def getValuesAtPoint(indir, rasterfileList, pos, lon, lat):
    #gt(2) and gt(4) coefficients are zero, and the gt(1) is pixel width, and gt(5) is pixel height.
    #The (gt(0),gt(3)) position is the top left corner of the top left pixel of the raster.
    for i, rs in enumerate(rasterfileList):

        presValues = []
        gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
        gt = gdata.GetGeoTransform()
        band = gdata.GetRasterBand(1)
        nodata = band.GetNoDataValue()

        x0, y0 , w , h = gt[0], gt[3], gt[1], gt[5]

        data = band.ReadAsArray().astype(numpy.float)
        #free memory
        del gdata

        if i == 0:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                Xc = x0 + x*w + w/2 #the cell center x
                y = int((p[1][lat] - y0)/h)
                Yc = y0 + y*h + h/2 #the cell center y
                try:
                    if data[y,x] != nodata:
                        presVAL = [p[1][lon],p[1][lat], '{:.6f}'.format(Xc), '{:.6f}'.format(Yc), data[y,x]]
                        presValues.append(presVAL)
                except:
                    pass
            df = pandas.DataFrame(presValues, columns=['x', 'y', 'Xc', 'Yc', rs])
        else:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                y = int((p[1][lat] - y0)/h)
                try:
                    if data[y,x] != nodata:
                        presValues.append(data[y,x])
                except:
                    pass
            df[rs] = pandas.Series(presValues)
    del data, band
    return df

Example of how to run it given that the rasters are in your current working directory:

rasDf = getValuesAtPoint('.', ['raster1', 'raster2'], inPoints, 'x', 'y')

Answer 10

If you have access to FME then you can use one of two transformers in FME Workbench.

The RasterCellCoercer ("Decomposes all input numeric raster features into individual points or polygons. One vector feature is output for each cell in the raster.")

The PointOnRasterValueExtractor ("Takes in point features and a single reference raster. The output consists of the band and palette value(s) at the location of each point.")

Answer 11

You can use Sample Raster Values from the Processing Toolbox in QGIS. This algorithm creates a new vector layer with the same attributes of the input layer and the raster values corresponding on the point location. If the raster layer has more than one band, all the band values are sampled.

Answer 12

Quick thought:

1. gdal_polygonize.py - polygonize your raster feature
2. Insert your point features and polygons into PostGIS database
3. Use st_intersects function to pull all elevation values where features intersect