Estimate cloudiness based on SCL layers at pixel level - Spatial Resampling needed

[VicenteYago]

At pixel level can be calculated yes, but due to the lesser resolution of the SCL layer (20mts/px) can not be compared at pixel level directly with the pixel values of the NDVI or SWIR values (10 mts/px)

[VicenteYago]

For the moment i will compute the global percentage of the image rather than pixel wise logical defective status:

def is_defective_px(SCL_px_val):

    if ((SCL_px_val == 4) | (SCL_px_val == 5)): # vegetated OR not_vegetated
        return False

    else : 
        return True

# Vectorized version
is_defective_px_vect = np.vectorize(is_defective_px)

# TODO: maybe in the denominator are included alpha values and thus
# computing an undervaluated percentage ?
# TODO type checks 
# TODO shape checks  
def S2_defective_px_perc(scl_raster) : 
    not_defective = sum(~np.array(is_defective_px_vect(scl_raster)))
    perc = not_defective / len(np.ndarray.flatten(scl_raster))
    return perc

As a result all the pixels of the global dataframe will have the same cloudiness_perc value

[VicenteYago]

The key problem is the 10m/px res bands such as :

• NDVI
• R
• NIR
• SWIR

have a shape of (1214, 2607) an thus when performing a 2x resampling from the SCL band (20mt/px):

upscale_factor = 2
# resample data to target shape
scl  = scl.read(1,
                out_shape=(
                    scl.count,
                    int(scl.height * upscale_factor),
                    int(scl.width * upscale_factor) 
                ), 
                resampling=Resampling.nearest,
                masked = True)

creates transformation a : (607, 1304) -> (1214, 2608), not matching the other layers exact height and width

Left: 20px/m SCL layer, Right: 10px/m resampled

[VicenteYago]

The issue if fixed performing a left join once the scl band has been supersampled to 10m/px with the gdal library:

from osgeo import gdal, gdalconst                   #

def resample_raster_gdal_nn(input_file, ref_file, out_file):
    # Opening input
    input = gdal.Open(input_file, gdalconst.GA_ReadOnly)
    inputProj = input.GetProjection()
    inputTrans = input.GetGeoTransform()

    # Opening ref
    reference = gdal.Open(ref_file, gdalconst.GA_ReadOnly)
    referenceProj = reference.GetProjection()
    referenceTrans = reference.GetGeoTransform()
    bandreference = reference.GetRasterBand(1)    
    x = reference.RasterXSize 
    y = reference.RasterYSize

    # Resampling
    driver= gdal.GetDriverByName('GTiff')
    output = driver.Create(out_file,x,y,1,bandreference.DataType)
    output.SetGeoTransform(referenceTrans)
    output.SetProjection(referenceProj)
    gdal.ReprojectImage(input,output,inputProj,referenceProj,gdalconst.GRA_NearestNeighbour)
    del output
    del input
    del reference

Usage :

scl_fp = './sen2r/out/SCL/S2A2A_20190120_141_Walnut-Gulch_SCL_10.tif'
ndvi_fp = './sen2r/indices/NDVI/S2A2A_20190120_141_Walnut-Gulch_NDVI_10.tif'
scl_10_fp = './sen2r/out/SCL_res10/20190120_SCL_10m_resampled_by_gdal.tif'
cd SCL_res10/
resample_raster_gdal_nn(input_file = scl_fp, ref_file = ndvi_fp, out_file = scl_10_fp)

In the add_scl_col function :

def add_scl_col(scl_fp, ndvi_fp, local_df, date,  scl_dir = "./sen2r/out/SCL_res10/"):

    ...

    scl_10_df = pd.DataFrame({
        'utm_x' : x,
        'utm_y' : y,
        'scl_value' : np.delete(scl_10_flatten, scl_10_flatten == 0)
    })
    scl_10_df = scl_10_df.astype('int')    
    local_df = pd.merge(local_df, scl_10_df, on = ["utm_x", "utm_y"], how = "left")
    return(local_df)

The left merge preserve the original pixels displaced in the resampled image, setting it with NaN. Hopefully the NaN values represent a negligible amount of the total, always affecting some of the pixels at the image margins.