Extract raster values at vector points location with Python/GDAL
по адресу http://gis-lab.info/qa/extract-values-rasters-eng.html
Script to extract raster values at points.
Extracting data from a single raster (or several rasters) using a set of vector point objects is a fairly common task. This is necessary when we want to construct the spectral profiles, quantify the accuracy of raster data, inspect changes or prepare training datasets for classification, etc.
Here is already a similar tool, on our website, but it was developed for ArcView and has several limitations.
Preparation
To use this tool you need our python script, Python 2.5 or higher and GDAL library with Python bindings.
The simplest way to get all necessary components under OS Windows is to use OSGeo4W installer:
- Download installer. If the previous link is not accessible, there is an alternative link on a mirror website
- Run the installer
- Select «Advanced install» and press «Next»
- Select gdal and gdal-python packages (under Libs category), then press «Next»
- Wait while requested packages and their dependencies are installed
If you run Linux, then install all necessary packages using your package manager. In Ubuntu/Debian you need to install python-gdal package.
Extract the script from archieve and copy extract_values.py into any directory in your PATH (for example, into C:\Tools\GIS or ~/bin). To run the script under Windows, you need to open OSGeo4W shell and run the script from the shell.
Usage
There is no graphical interface for the script and it runs in console mode at the moment. It takes several arguments as input and output
extract_values.py [-c] [-r] point_shapefile [raster_file(s)] [-d directory_with_rasters]
- point_shapefile — vector point in ESRI shapefile format. Raster values will be extracted at these locations
- raster_file(s) — raster(s) from which values will be extracted. If rasters and script located in different directories, full path should be specified. Use space to separate multiple rasters
- -c — write data into external file. A CSV file (with same name as shapefile) will be created within the input shapefile folder, in this file all attributes and extracted data will be written. In this mode input shapefile remains unchanged
- -r — use this switch to reproject points into raster CRS, when rasters and vector are in different CRS
- -d — directory processing mode. After this switch a path to directory with rasters should be specified. Data will be extarcted from all rasters in this directory
The script works with single- and multiband raster in all GDAL-compatible formats. For each single-band raster in point shapefile, a new field will be created. The name of this field is is identical to raster name. For multiband rasters the field name in DBF file will consists of raster name and band number.
Important! You must remember that the field name in the DBF format is limited to 10 characters, hence long names will be truncated.
Here is result of the data extraction from 6-bands Landsat fragment (raster data type is Byte, raster filename — clearcuts_174016.tiff), opened in QGIS. For the initial two fields (id, taxon) new fields clearcut_1 - clearcut_6 (see, that filed names are truncated due to the DBF format limitations) with corresponding values were added.\
Examples
Extract values from raster after.tiff using points from shapefile poi.shp (script and files are in the same directory)
extract_values.py poi.shp after.tiff
Extract values from raster after.tiff using points from shapefile poi2.shp and write result into external file
extract_values.py -с poi2.shp after.tiff
Extract values from raster after.tiff using points from shapefile poi2.shp, raster and vector are in different CRS (so we use -r switch to enable reprojection)
extract_values.py -r poi2.shp after.tiff
Extract values from rasters before.tiff and after.tiff using points from shapefile points.shp
extract_values.py points.shp before.tiff after.tiff
Extact values from all rasters in directory using points from shapefile points.shp
extract_values.py points.shp -d D:\data\rasters_veg