11

I'm a beginner with python. I would like to extract raster (.tif) values from a set of specific points that I have in a shapefile. I would like to have the same result as the QGIS Point Sampling Tool plugin. I can't install on Python libraries like qgis.core arcpy pygrass, but I have available GeoPandas GDAL Fiona Shapely rasterio and others.

How can I do this?

3
  • 1
    Did you look at rasterstats?
    – Jon
    Apr 1, 2019 at 17:15
  • 1
    Have a look here : gis.stackexchange.com/a/310304/93097 Apr 1, 2019 at 17:27
  • 1
    Hum... You should try simple things first, like installing QGIS. As @ndawson pointed out, QGIS can generate the code for you.
    – jgrocha
    Jun 4, 2019 at 23:09

3 Answers 3

23

I use rasterio and geopandas. My example uses UTM coordinates. Obviously these fields will depend on your particular shapefile. In my experience this produces indentical results to the QGIS Point Sampling Tool. I like this method because the resulting DataFrame of point and corresponding raster values is easy to analyze (e.g. compute the difference between the point and raster values) and then plot or export to some other tabular format (e.g. CSV).

import rasterio
import geopandas as gpd

# Read points from shapefile
pts = gpd.read_file('your_point_shapefile.shp')
pts = pts[['UTM_E', 'UTM_N', 'Value', 'geometry']]
pts.index = range(len(pts))
coords = [(x,y) for x, y in zip(pts.UTM_E, pts.UTM_N)]

# Open the raster and store metadata
src = rasterio.open('your_raster.tif')

# Sample the raster at every point location and store values in DataFrame
pts['Raster Value'] = [x for x in src.sample(coords)]
pts['Raster Value'] = probes.apply(lambda x: x['Raster Value'][0], axis=1)
7
  • 11
    You can avoid the last line if you do pts['Raster Value'] = [x[0] for x in src.sample(coords)] in the line above. Dec 16, 2019 at 11:20
  • 2
    If you don't have specific coordinates stored in your shapefile's fields, or just want more generic code, you can instead use: zip(pts.geometry.x, pts.geometry.y) (apologies for the lack of formatting, none of the instructions on formatting actually work when I use them)
    – Anders
    Mar 6, 2020 at 15:28
  • 2
    thats awesome!! But what exactly is this probes?
    – Lenn
    Oct 17, 2020 at 10:09
  • @Lenn Probes is the variable name for the point locations where I same the raster. In my case they were snow depth probes. Oct 19, 2020 at 17:46
  • What interpolation is used? Or is it a nearest neighbor sampling? Mar 9 at 1:37
1

You can use the "Sample raster values" tool from the Processing toolbox. Like any of the Processing tools, if you run the tool from the GUI, and then look at the Processing history window you can then click the corresponding entry in the log, you can get an equivalent PyQGIS command which does the same operation:

enter image description here

1
  • 1
    Like I said before I have problem to install Qgis's library, so I can't use this type of algorithm
    – kino
    Apr 2, 2019 at 7:55
0

I am sure that it can be done in a more pythonic way, but I use a simple python script, which calls gdallocationinfo process and exports results to CSV (you must have GDAL installed). If you want to try, just copy the script below, fill in your sites list and path to your raster file.

#!/usr/bin/env python
#created by jurajb
import os
from subprocess import Popen,PIPE
import csv

sites= [
['Point1',-15.39495,28.33711],
['Point2',-15.548307,28.248216]
]  ## 'Name',Lat,Lon

rast = '/path/to/raster/my_raster.tif'
param = 'Raster_VALUE'
csvoutfl = 'sites.csv'

##==

scr = open(csvoutfl, 'w')
header = 'site,lat,lon,{}'.format(",".join(param))

for i in sites:
    csvline = '{},{},{}'.format(i[0],i[1],i[2])
    result = os.popen("gdallocationinfo -wgs84 -valonly {0} {1} {2}".format(rast, i[2], i[1])).read()
    try:
        result = float(result)
    except ValueError:
        result = 'Err'
    csvline += ',{}'.format(result)
    scr.write('{}\n'.format(csvline))
scr.close()
print "\n\nCREATED:  {}\n\n=== Finished ===".format(csvoutfl)

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.