I am trying to write a script to perform some zonal statistics with rasterestats, but I am really struggling with nodata values. I noticed a different behaviour 1) when running zonal_stats
with numpy array as raster input (doesn't ignore nodata) and 2) when I specify the path to the raster directly (it works and ignores nodata). But I need to include the process in a longer script which requires to use the first option.
I have tried both reading rasterio Dataset with and without masked=True
option, setting various values for nodata in the original raster (including: None, 0, -3.4028234663852886e+38) and using numpy ma.masked_where/ma.masked_values
on the array but the result is always the same, with the corresponding nodata value showing up in the output statistics. Any explanation and workaround for this?
import numpy as np
import time, rasterio
from rasterio.plot import show
from rasterstats import zonal_stats
import geopandas as gpd
# crearte function to print out raster metadata
def rasterinfo(r):
print('ObjectType: ',type(r))
kwds=r.profile
print('Driver: ',r.driver)
print('CRS: ',r.crs)
print('Width: ',r.width)
print('Height: ',r.height)
print('Number of bands: ',r.count)
print('Data type:',kwds['dtype'])
print('Boundaries: ',r.bounds)
print('No data values: ',kwds['nodata'])
print('Spatial.Res:',int(r.transform.a))
print('Summary Profile: \n',r.profile)
# set path to raster and zone vector
rst = r'F:\GeoData\raster\clc_resample_02.tif'
fp= r'F:\GeoData\Vector\nuts3_2013_suppid.shp'
# read in raster with rasterio and check metadata
src = rasterio.open(rst)
rasterinfo(src)
# read zone vector to geodataframe and set crs
zone = gpd.read_file(fp)
zone = zone[['nutsuppid','geometry']].to_crs(src.crs)
arr = src.read(1, masked=True)
affine= src.transform
# perform zonal statistics
zs = zonal_stats(zone, arr, affine=affine)[0:3] # filtering for testing purposes
zs
Output:
ObjectType: <class 'rasterio.io.DatasetReader'>
Driver: GTiff
CRS: EPSG:3035
Width: 6500
Height: 4600
Number of bands: 1
Data type: int16
Boundaries: BoundingBox(left=900000.0, bottom=899999.9999999991, right=7400000.000000001, top=5500000.0)
No data values: -32768.0
Spatial.Res: 1000
Summary Profile:
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32768.0, 'width': 6500, 'height': 4600, 'count': 1, 'crs': CRS.from_epsg(3035), 'transform': Affine(1000.0000000000001, 0.0, 900000.0,
0.0, -1000.0000000000001, 5500000.0), 'tiled': False, 'interleave': 'band'}
[{'min': -32768.0, 'max': 240.0, 'mean': -17119.210227272728, 'count': 704},
{'min': -32768.0, 'max': 240.0, 'mean': -12155.484069312464, 'count': 1789},
{'min': -32768.0, 'max': 240.0, 'mean': -15936.236698499319, 'count': 1466}]
test = zonal_stats(zone, rst)[0:3]
test
Output:
[{'min': 212.0, 'max': 240.0, 'mean': 216.2754491017964, 'count': 334},
{'min': 212.0, 'max': 240.0, 'mean': 215.7128801431127, 'count': 1118},
{'min': 212.0, 'max': 240.0, 'mean': 220.4558823529412, 'count': 748}]
I check the masked array
arr
masked_array(
data=[[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
...,
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --],
[--, --, --, ..., --, --, --]],
mask=[[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True]],
fill_value=-32768,
dtype=int16)
Finally, if I add the categorical=True
option to the zonal statistics function I get similar results, confirming that the process somehow works in both cases but it doesn't ignore nodata when using the array (my testing raster has categorical data between 212 and 240)
zs = zonal_stats(zone, arr, affine=affine, categorical=True)[0:3]
zs
Output:
[{-32768.0: 370, 212.0: 264, 221.0: 28, 240.0: 42},
{-32768.0: 671, 212.0: 851, 221.0: 175, 240.0: 92},
{-32768.0: 718, 212.0: 516, 221.0: 9, 240.0: 223}]
test = zonal_stats(zone, rst, categorical=True)[0:3]
test
Output:
[{212: 264, 221: 28, 240: 42},
{212: 851, 221: 175, 240: 92},
{212: 516, 221: 9, 240: 223}]