20

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.

  • 1
    So you have points and you need to extract the values from the raster under those points, or do you need to convert the raster cells into points. Just checking before I try and work out the answer. – Nathan W Nov 16 '10 at 11:07
  • The first, i have the points and i need to extract thevalues from the raster, under those points. THNX!! – Vassilis Nov 16 '10 at 12:50

11 Answers 11

35

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/

  • 5
    People are still finding the above-mentioned post through this Q&A. However, the plugin is not the only solution, and not always the easiest solution anymore. An alternative solution is the Saga function ‘Add grid values to point’ in the processing toolbox. See for details this post. – Ecodiv Jul 4 '14 at 10:17
  • Caution. I just ran the Point Sampling Tool. 60,000 points and 13 rasters. The results failed my 30 random sample test for each year. This tool has issues with large datasets. I would not use it. – If you do not know- just GIS Feb 11 '16 at 15:29
  • Despite the mentioned issues with large datasets, it is very useful for extracting all multiband values in one go. All other QGIS-related solutions do only support the extraction of one band (like GRASS r.what) or prohibit the use of multiband raster (like Saga - raster values to points). – EikeMike Jun 6 '17 at 12:53
6

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).

4

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

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.

3

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

  • 2
    For Grass I found this:pvanb.wordpress.com/2010/05/05/… – Vassilis Nov 19 '10 at 8:58
  • Another GRASS module that is capable of sampling a raster is v.sample, which is also available through the QGIS processing toolbox. – user55937 Feb 12 '16 at 20:45
3

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.

3

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))
  1. 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.

  • 1
    this is the best answer for this question – B-C B. Jun 25 at 15:38
2

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 :)

1

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.")

  • No I don't have or use FME, is a standalone application or a plugin? – Vassilis Nov 26 '10 at 9:14
1

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')
0

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
  • interesting approach, because yesterday start studing how to use Postgis. – Vassilis Mar 12 '12 at 20:35
  • Thanks, it's fairly simplistic but it works. Here is what I was able to produce with this approach: i.imgur.com/h8CGJ.png – user5584 Mar 14 '12 at 13:58

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.