Extract point value from a raster file with Python, Geopandas and Rasterio - Tutorial
/This tutorial has a complete case of spatial analysis for the extraction of point data from a raster dataset with Python and its libraries Geopandas and Rasterio. The procedure is entirely geospatial and uses shapefiles and tifs as input data; data calculation was performed on a Jupyter Lab environment.
Tutorial
Script
This is the complete Python script:
#import required libraries
%matplotlib inline
import matplotlib.pyplot as plt
import geopandas as gpd
import rasterio
from rasterio.plot import show
C:\Users\GIDA2\Anaconda3\lib\site-packages\geopandas\_compat.py:88: UserWarning: The Shapely GEOS version (3.4.3-CAPI-1.8.3 r4285) is incompatible with the GEOS version PyGEOS was compiled with (3.8.1-CAPI-1.13.3). Conversions between both will be slow.
shapely_geos_version, geos_capi_version_string
#open point shapefile
pointData = gpd.read_file('Shp/pointData.shp')
print(pointData.crs)
pointData.plot()
epsg:32611
<matplotlib.axes._subplots.AxesSubplot at 0x28cb7d97f08>
#open raster file
ndviRaster = rasterio.open('Rst/ndviImage.tiff')
print(ndviRaster.crs)
print(ndviRaster.count)
EPSG:32611
1
#show point and raster on a matplotlib plot
fig, ax = plt.subplots(figsize=(12,12))
pointData.plot(ax=ax, color='orangered')
show(ndviRaster, ax=ax)
<matplotlib.axes._subplots.AxesSubplot at 0x28cb90cf248>
#extract xy from point geometry
for point in pointData['geometry']:
print(point.xy[0][0],point.xy[1][0])
239728.82443269365 3990037.7658296037
249272.44019424845 3995708.972054402
247098.89419619803 3987464.4872342106
254893.6798443789 3991661.679506308
255893.01133773543 3997207.9692944367
257491.94172710588 3984091.7434441326
243001.63507343628 3978745.319954675
238054.94418132148 3978370.5706446664
#extract point value from raster
for point in pointData['geometry']:
x = point.xy[0][0]
y = point.xy[1][0]
row, col = ndviRaster.index(x,y)
print("Point correspond to row, col: %d, %d"%(row,col))
print("Raster value on point %.2f \n"%ndviRaster.read(1)[row,col])
Point correspond to row, col: 708, 332
Raster value on point 0.11
Point correspond to row, col: 519, 650
Raster value on point 0.34
Point correspond to row, col: 794, 578
Raster value on point 0.47
Point correspond to row, col: 654, 837
Raster value on point 0.26
Point correspond to row, col: 469, 871
Raster value on point 0.33
Point correspond to row, col: 906, 924
Raster value on point 0.44
Point correspond to row, col: 1084, 441
Raster value on point 0.26
Point correspond to row, col: 1097, 276
Raster value on point 0.36
Input data
You can download the input data from this link.