Extract raster values to CSV in Python

Question

Attempting to gather MODIS aerosol data for a number of different AERONET stations in the UK to plot against one another. However, due to the number of files (upwards of 2000), QGIS and ArcMap crash.

I have attempted to find a solution within RSGISLib, whereby I extract a raster to .csv via a defined vector but keep getting the following:

ERROR 1: Attempt to read shape with feature id (1) out of available range. Segmentation fault

I have seen a method to extract to a .csv using r, however, I would prefer for this to be done within Python. Is there a way?

import rsgislib
from rsgislib import zonalstats
import glob

inputimage = '/home/student/Desktop/Ben_Folder/Test2/test.kea'
inputvector = '/home/student/Desktop/Ben_Folder/AERONET_POINTS/Chilbolton.shp'
outputtxt = 'test.csv'
useBandNames = False
zonalstats.pointValue2TXT(inputimage, inputvector, outputtxt, useBandNames)

Answer 1

I would do it with Rasterio and Fiona. In the following code snippet, I use pandas to create a Serie object and write the records to a CSV file in order to improve code readability ; however you could also use the csv module. Also, I am assuming that your shapefile contains points and not polygons.

import pandas as pd
import rasterio
import fiona

values = pd.Series()

# Read input shapefile with fiona and iterate over each feature
with fiona.open('stations.shp') as shp:
    for feature in shp:
        station_name = feature['properties']['name']
        coords = feature['geometry']['coordinates']
        # Read pixel value at the given coordinates using Rasterio
        # NB: `sample()` returns an iterable of ndarrays.
        with rasterio.open('test.kea') as src:
            value = [v for v in src.sample([coords])][0][0]
        # Update the pandas serie accordingly
        values.loc[station_name] = value

# Write records into a CSV file
values.to_csv('test.csv')