Numpy is not setting properly nan values: arr[arr== 0] = np.nan

Question

I am creating a script to normalize a satellite scene. My input image is of type float32, and no NoData value is assigned. However, during the normalization, I want to avoid using pixels with a value of 0 (usual black borders in the scene). To do that, I do the following for each band: current_band[current_band == 0] = np.nan, but this assigns nan values to all pixels instead of only those with 0.

I want to avoid setting the NoData value = 0 to the input file to avoid the output file inheriting that property. And even if I set it, numpy it is not considering it.

Any idea of what am I doing wrong?

Here you are my code:

import numpy as np
import rasterio
import argparse

# TODO: implement check if file is float32

def normBands(input_raster, output_raster):
    with rasterio.open(input_raster) as src:
        # Read the number of bands and the dimensions of the image
        num_bands = src.count
        height = src.height
        width = src.width
        print(src.nodata) # this should print None

        # Create an empty array to store the normalized bands
        norm_bands = np.empty((num_bands, height, width), dtype=np.float32)
        norm_bands[:] = np.nan

        # Loop over all the bands
        for band in range(num_bands):
            # Read the current band
            current_band = src.read(band + 1)
            current_band[current_band == 0] = np.nan #NOT WORKING, makes all pixels as nan

            # Calculate mean and std for the current band
            mean = np.mean(current_band)
            std = np.std(current_band)
            print("band: "+ str(band+1))
            print(mean) ## smth NOT WORKING, prints nan!!!
            print(std)  ## smth NOT WORKING, prints nan!!!

            # Normalize the current band mean and std
            norm_band = (current_band - mean)/std

            # Add the normalized band to the array of normalized bands
            norm_bands[band, :, :] = norm_band

        # Save the normalized image
        with rasterio.open(output_raster, 'w', **src.profile) as dst:
            dst.write(norm_bands)

if __name__ == "__main__":
    aparser = argparse.ArgumentParser(
        description='normalize bands'
    )
    aparser.add_argument(
        '--input-raster',
        default=None
    )
    aparser.add_argument(
        '--output-raster',
        default=None
    )
    args = aparser.parse_args()
    normBands(args.input_raster, args.output_raster)

Answer 1

Your arrays are not all nan.

It prints nan because that is how np.mean and np.std work, if the array contains any nans, the result will be nan.

You can use nanmean and nanstd instead:

import numpy as np
a = np.array([[1, 2, 3],
              [4, 5, np.nan],
              [np.nan, 6, np.nan],
              [np.nan, np.nan, np.nan]])

print(np.mean(a))
#nan

print(np.nanmean(a))
#3.5

print(np.std(a))
#nan

print(np.nanstd(a))
#1.707825127659933