Multiply arrays from two directories to get cloudfree multiband raster

JohannesWissing — Sun, 12 Dec 2021 00:58:58 GMT

I would like to multiply 4band raster from a directory with corresponding 1 band cloud masks to get cloudfree 4band raster in the end. In my attempt, I started with multiplying the first band by the cloud mask. The cloud mask raster are already reclassified ( pixels with values = 0, pixels with no values or clouds= nodata). The problem is that the scripts returns just sometimes one file and returns errors I don't understand:

ERROR 1: Deleting D:/SR\20170207_4band.tif failed:Permission denied

-----------

import pandas as pd
import gdal
import os
import numpy as np
import glob

SR_images = 'D:/SR'
sr = glob.glob(os.path.join(SR_images, '*.tif'))
reclass = 'D:/reclass'
mask = glob.glob(os.path.join(reclass, '*.tif'))
savedest = r'D:\maskedsr'

for e, i in zip(sr, mask):
   srraster = gdal.Open(os.path.join(SR_images, e))
   maskraster = gdal.Open(os.path.join(reclass, i))
   band1 = srraster.GetRasterBand(1).ReadAsArray()
   band2 = srraster.GetRasterBand(2).ReadAsArray()
   band3 = srraster.GetRasterBand(3).ReadAsArray()
   band4 = srraster.GetRasterBand(4).ReadAsArray()
   bandmask = maskraster.GetRasterBand(1).ReadAsArray()
   geotransform = srraster.GetGeoTransform()
   geoproj = srraster.GetProjection()
   xsize = srraster.RasterXSize
   ysize = srraster.RasterYSize
   format1 = "GTiff"
   driver = gdal.GetDriverByName(format1)
   dst_datatype = gdal.GDT_UInt16

   cal_band1 = band1 * bandmask
   name = e
   output_path = os.path.join(savedest, name)
   dst_ds = driver.Create(output_path, xsize, ysize, 1, dst_datatype)
   dst_ds.GetRasterBand(1).WriteArray(cal_band1)
   dst_ds.SetGeoTransform(geotransform)
   dst_ds.SetProjection(geoproj)
   dst_ds.FlushCache()
   dst_ds = None

Thanks!

topic Multiply arrays from two directories to get cloudfree multiband raster in Python Questions

Multiply arrays from two directories to get cloudfree multiband raster