centreborelli / s2p

Satellite Stereo Pipeline
GNU Affero General Public License v3.0
208 stars 67 forks source link

GML ROI mask not working properly #117

Open daviddemeij opened 2 years ago

daviddemeij commented 2 years ago

For some reason when I use an ROI GML file as roi parameter the tiles don't get properly masked. I believe the problem occurs in the cldmask code. When I rasterize the GML file and input it as wat parameter the tiles are masked properly.

When I read the mask.png files for each tile when using the GML file as roi input all the values are 1.

Or maybe I am doing something wrong? I am using Pleiades images and setting the roi to IMG_PHR1A_P_001/MASKS/ROI_PHR1A_P_202106151554566_SEN_5789776101-1_MSK.GML.

example Result when using the GML ROI mask:

Result when using the rasterized ROI mask:

A potential solution could be to always first rasterize any GML input file. I used the following code to rasterize the GML file

import rasterio
import rasterio.features
import numpy as np
import geopandas as gpd

def rasterize_gml_mask(roi_gml_path, raster_path, out_path, cld_gml_path=None):
    """
    Rasterize a GML mask into a raster.

    Args:
        roi_gml_path (str): Path to the GML file containing the ROI.
        raster_path (str): Path to the raster to rasterize the ROI into.
        out_path (str): Path to the output raster.
        cld_gml_path (str): Path to the GML file containing the cloud mask.
    """
    # Set ROI pixels to 1
    shapes = [(s, 1) for s in gpd.read_file(roi_gml_path).geometry.tolist()]
    if cld_gml_path:
        try: 
            cloud_shapes = gpd.read_file(cld_gml, driver='gml').geometry.tolist()
        except ValueError as e:
            if str(e) == "Null layer: ''":
                print(f"No clouds, skipping rasterizing clouds")
                cloud_shapes = []
            else:
                raise e
        # Set cloud pixels to 0
        shapes += [(s, 0) for s in cloud_shapes]

    with rasterio.open(raster_path) as ds:
        out = np.zeros(ds.shape, dtype=np.uint16)

        # Rasterize the shapes
        rasterio.features.rasterize(shapes, out=out)

        # Write rasterized mask to out_path
        with rasterio.open(out_path, 'w', driver='GTiff', count=1, dtype=np.uint16, width=ds.width, height=ds.height) as dst:
            dst.write(out, 1)

However this changes the requirements of s2p pipelines