I wrote a code based from rasterio and numpy that gets the nanmean of every pixel from multiple of raster images. I first read it using rasterio, then I convert the array to float32 so that I can set the NaN (-32768) values to np.nan. I then add each image to a list.

When the list is filled. I then compute it using np.nanmean from numpy. Real enough that NaN values were ignored in computing the output image. However, I do not understand the values from the resulting image. The input raster images have a value from 0 to 5. My output has a value 0 - 20000.

Note: The input images are really Float32 in QGIS. I dont know why it was converted to Int16 in rasterio.

  1. Why was it converted to Int16?
  2. How do I correct this? I want an output with nanmean from the inputs. The real values should be between 0 - 5 also.
import rasterio
from glob import glob
import numpy as np

path = '/my/directory/*.nc'
output = 'output/AOT_L2_Mean_MAM.tif'
datasets = glob(path)

rasters = []
for data in datasets:
    ds = rasterio.open(f'netcdf:{data}:AOT_L2_Mean')
    band = ds.read(1)
    band = band.astype('float32')
    band[band==-32768] = np.nan
out = np.nanmean(rasters, axis=0)

with rasterio.open(
    crs='EPSG:4326', # +proj=latlong
) as dst:
    dst.write(out, 1)

1 Answer 1


Have you tried:

band = ds.read(1, masked=True)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.