sentinel-hub / sentinel2-cloud-detector

Sentinel Hub Cloud Detector for Sentinel-2 images in Python
Creative Commons Attribution Share Alike 4.0 International
428 stars 93 forks source link

[HELP] Cloud detector run locally is masking non-clouds #80

Closed eoasker closed 9 months ago

eoasker commented 9 months ago

Hello. I am trying to use the Cloud-detector on locally stored S2-Images in .SAFE format (L1C). I ran the cloud detector for 3 images. One of them masked quite well the clouds. However when I ran the same code for two other images with almost no clouds in the image, the cloud mask that is returned comes completely wrong. Is there something wrong in my code? Or is this a problem with the detection?

This is the code I am running:

        #run s2cloudless
        #initialize s2cloudless
        cloud_detector = S2PixelCloudDetector(threshold=0.4, average_over=4, dilation_size=2)
        print("Resampling bands to 10 meters...")
        with rasterio.open(f"{output_clip}/{jp2_name}_B02.jp2") as img:
            B02 = img.read()
            tmparr = np.empty_like(B02)
            aff = img.transform
        with rasterio.open(f"{output_clip}/{jp2_name}_B04.jp2") as img:
            B04 = img.read()         
        with rasterio.open(f"{output_clip}/{jp2_name}_B08.jp2") as img:
            B08 = img.read()
        with rasterio.open(f"{output_clip}/{jp2_name}_B01.jp2") as img:
            B01 = img.read()
            reproject(
                B01, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B01 = tmparr            
        with rasterio.open(f"{output_clip}/{jp2_name}_B05.jp2") as img:
            B05 = img.read()
            reproject(
                B05, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B05 = tmparr            
        with rasterio.open(f"{output_clip}/{jp2_name}_B8A.jp2") as img:
            B8A = img.read()
            reproject(
                B8A, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B8A = tmparr            
        with rasterio.open(f"{output_clip}/{jp2_name}_B09.jp2") as img:
            B09 = img.read()
            reproject(
                B09, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B09 = tmparr              
        with rasterio.open(f"{output_clip}/{jp2_name}_B10.jp2") as img:
            B10 = img.read()
            reproject(
                B10, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B10 = tmparr            
        with rasterio.open(f"{output_clip}/{jp2_name}_B11.jp2") as img:
            B11 = img.read()
            reproject(
                B11, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B11 = tmparr             
        with rasterio.open(f"{output_clip}/{jp2_name}_B12.jp2") as img:
            B12 = img.read()
            reproject(
                B12, tmparr,
                src_transform = img.transform,
                dst_transform = aff,
                src_crs = img.crs,
                dst_crs = img.crs,
                resampling = Resampling.bilinear)
            B12 = tmparr                   

        bands = np.array([np.dstack((B01[0]/10000.0,B02[0]/10000.0,B04[0]/10000.0,B05[0]/10000.0,B08[0]/10000.0,B8A[0]/10000.0,B09[0]/10000.0,B10[0]/10000.0,B11[0]/10000.0,B12[0]/10000.0))])

        #run s2cloudless
        print("Running s2cloudless..")
        cloud_probs = cloud_detector.get_cloud_probability_maps(bands)
        mask = cloud_detector.get_cloud_masks(bands).astype(rasterio.uint8)

        print("Saving cloud mask as raster..")
        #save clouds mask
        with rasterio.open(f"{output_mask}/cloud_mask.tif", "w",  driver='GTiff',height=mask.shape[1],width=mask.shape[2],count=1,dtype=rasterio.uint8,crs=CRS,transform=aff) as dest:
            dest.write(mask)
        with rasterio.open(f"{output_mask}/cloud_probability_mask.tif", "w",  driver='GTiff',height=cloud_probs.shape[1],width=cloud_probs.shape[2],count=1,dtype=cloud_probs.dtype,crs=CRS,transform=aff) as dest:
            dest.write(cloud_probs)

The image I am showing you below is _S2A_MSIL1C_20231004T113321_N0509_R080_T29TNF20231004T151316.SAFE and I have clipped it around the district of Porto. This is the RGB composite and the mask cloud for the same image: image

image

mlubej commented 9 months ago

Hi @eoasker,

can you provide some info, please:

Are you plotting masks or cloud probabilities above a certain threshold?

eoasker commented 9 months ago

Hi @mlubej . So the bands shape is (1, 5210, 7660, 10) and the mask shape is (1, 5210, 7660) . I am opening the mask and cloud probability on qgis, where for the mask the min and max is 0 and 1 which is the image I am showing you above.

mlubej commented 9 months ago

OK, it looks like the shapes are correct. Can you use the info tool on the black and white areas to check the values of the tiff if they are correct? I'm suspicious of this being a no_data issue, where the 0 values of the mask were tagged as no_data automatically.

Try setting the nodata to -1 in qgis manually, or add a specific value when you save the tiff.

Let me know if this solves the issue.

eoasker commented 9 months ago

I have solved the problem. While I was dividing the bands by 10000.0 I also needed to remove the offset from them as seen here . I also made sure that there were no negative values. Thank you for the help.