Exporting multiband Earthpy masked GeoTIFF with Rasterio resulted in original non-masked file

Question

I currently wanted to mask a stacked Landsat 8 GeoTIFF file with the Pixel QA band to remove the clouds and cloud shadows. So far, I have successfully followed the tutorial here, and have properly plotted the masked scene using a combination of EarthPy and Rasterio. The resulting masked scene is a NumPy masked array.

However when I try to export the clean array as a GeoTIFF file (clean.tif), the resulting file contains the original scene instead of the masked scene file. Below is the code that I have:

from rasterio.plot import plotting_extent
import rasterio as rio
import earthpy.plot as ep
import earthpy.mask as em

# Open mask file
with rio.open("mask.tif") as mask_src:
  mask = mask_src.read()
  mask_meta = mask_src.profile

# Open scene file
with rio.open("scene.tif") as scene_src:
  scene = scene_src.read()
  scene_ext = plotting_extent(scene_src)
  scene_trf = scene_src.transform
  scene_meta = scene_src.profile

# Perform masking
clean = em.mask_pixels(scene, mask)

# Print metadata for mask and scene files
print("masked scene shape => " + str(clean.shape))
print("mask meta => " + str(mask_meta))
print("scene meta => " + str(scene_meta))

# Open and write destination tif file
with rio.open("clean.tif", 'w', **scene_meta) as clean_dst:
  clean_dst.write(clean)

# Open destination tif file
with rio.open("clean.tif") as final_dst:
  final = final_dst.read()
  final_ext = plotting_extent(scene_src)

# Plot mask file
ep.plot_bands(mask)
# Plot scene file
ep.plot_rgb(scene, rgb=[4, 3, 2], extent=scene_ext, stretch=True)
# Plot masked scene
ep.plot_rgb(clean, rgb=[4, 3, 2], extent=scene_ext)
# Plot destination tif file
ep.plot_rgb(final, rgb=[4, 3, 2], extent=final_ext, stretch=True)

And here are the plots of the files:

1. mask file:
   

2. scene file:
   

3. plotted masked array:
   

4. saved masked array to geotiff:
   

Here is the Dropbox link to the files and the script. Really scratched my head around this, I appreciate any pointers as to what happened and how to fix it. Thank you all:)

Answer 1

Fixed it, I need to supply the NoData value that is taken from the original scene to the masked array. Here is the fixed script:

from rasterio.plot import plotting_extent
import rasterio as rio
import earthpy.plot as ep
import earthpy.mask as em

# Open mask file
with rio.open("mask.tif") as mask_src:
  mask = mask_src.read()
  mask_meta = mask_src.profile

# Open scene file
with rio.open("scene.tif") as scene_src:
  scene = scene_src.read()
  scene_ext = plotting_extent(scene_src)
  scene_trf = scene_src.transform
  scene_meta = scene_src.profile
  scene_nodata = scene_src.nodata

# Perform masking
clean = em.mask_pixels(scene, mask)
clean = clean.filled(scene_nodata)

# Print metadata for mask and scene files
print("masked scene shape => " + str(clean.shape))
print("mask meta => " + str(mask_meta))
print("scene meta => " + str(scene_meta))

# Open and write destination tif file
with rio.open("clean.tif", 'w', **scene_meta) as clean_dst:
  clean_dst.write(clean)

# Open destination tif file
with rio.open("clean.tif") as final_dst:
  final = final_dst.read()
  final_ext = plotting_extent(scene_src)

# Plot mask file
ep.plot_bands(mask)
# Plot scene file
ep.plot_rgb(scene, rgb=[4, 3, 2], extent=scene_ext, stretch=True)
# Plot masked scene
ep.plot_rgb(clean, rgb=[4, 3, 2], extent=scene_ext)
# Plot destination tif file
ep.plot_rgb(final, rgb=[4, 3, 2], extent=final_ext, stretch=True)

Hopefully this will be useful for people who are facing the same issue.