rasterstats zonal statistics does not ignore nodata

Question

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}]

Answer 1

Solved! I realized that I was changing only the meta property not the actual nodata value. I still don't understand why the masked array did not work but using the code below did the trick, just I had to convert to float in order to use np.nan. There might also be better solutions not requiring the conversion.

with rasterio.open('example.tif') as src:
    ndval = src.nodatavals[0]
    array = src.read(1)
    array = array.astype('float64')
    array[array==ndval] = np.nan