get_cloud_masks results into raster

[minarikro]

How can I save the result of get_cloud_masks binary numpy array into geotiff to work with mask layer in GIS?

[lvhengani]

@minarikro , I am assuming you are using cloudless with downloaded data (without using Sentinel Hub) like I do and have gdal installed.

from osgeo import gdal_array

cloud_detector = S2PixelCloudDetector(threshold=0.4, average_over=4, dilation_size=2)
cloud_probs = cloud_detector.get_cloud_probability_maps(my_band_images)
cloud_masks = cloud_detector.get_cloud_masks(my_band_images)

template_file = "/path/to/scene.SAFE/GRANULE/{granule}/IMG_DATA/band.jp2"

output = gdal_array.SaveArray(cloud_probs[0, :, :], out_file_prob, format="GTiff", prototype=template_file) 
output = gdal_array.SaveArray(cloud_masks[0, :, :], out_file_mask, format="GTiff", prototype=template_file) 

In my case the template_file is a geotif, I am not sure how it behaves with a jp2 file.

[minarikro]

@lvhengani thank you. I used sentinel hub, but I found another post describing it. I modified it for my purpose. But now I want to run s2cloudless to single downloaded image with all bands (in SAFE format?). Can you share with me some code example of your workflow? I did something wrong , because when I write outtput mask with rasterio, everything is 1.

Here is the code to write mask as geotiff with rasterio from sentinel hub wms request.

import rasterio
wms_bands_request = WmsRequest()
wms_bands = wms_bands_request.get_data()

cloud_detector = S2PixelCloudDetector(threshold=0.2, all_bands=True, average_over=4, dilation_size=2)
cloud_prob = cloud_detector.get_cloud_probability_maps(np.array(wms_bands))
cloud_mask = cloud_detector.get_cloud_masks(np.array(wms_bands)

dst_shape = cloud_mask[0]
dst_transform = rasterio.transform.from_bounds(*bounding_box.get_lower_left(), *bounding_box.get_upper_right(), width=dst_shape.shape[1], height=dst_shape.shape[0])

with rasterio.open("path", 'w', driver='GTiff', width=dst_shape.shape[1], height=dst_shape.shape[0], count=1, dtype=np.uint8, nodata=255, transform=dst_transform, crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") as dst: dst.write(dst_shape.astype(np.uint8), indexes=1)

[AleksMat]

Closing this for now. Feel free to reopen any time.