How to get nanmean of each pixel in multiple raster image?

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)
    #print(band.dtype)
    band = band.astype('float32')
    band[band==-32768] = np.nan
    rasters.append(band)
    
out = np.nanmean(rasters, axis=0)

with rasterio.open(
    output,
    'w',
    driver='GTiff',
    height=out.shape[0],
    width=out.shape[1],
    count=1,
    dtype=out.dtype,
    crs='EPSG:4326', # +proj=latlong
) as dst:
    dst.write(out, 1)

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