Skip to content
Snippets Groups Projects
Commit 087e39be authored by pascal.mouquet_ird.fr's avatar pascal.mouquet_ird.fr
Browse files

OTB removed from cloud masks processing

parent 5b988340
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
requirements.txt
+ 1
− 0
View file @ 087e39be
......@@ -11,4 +11,5 @@ typing>=3.6.6
pandas>=0.23.4
urllib3>=1.22
geopandas>=0.5.0
scipy
sen2chain/cloud_mask.py
+ 165
− 87
View file @ 087e39be
......@@ -8,13 +8,14 @@ import logging
import pathlib
from pathlib import Path
import pyproj
import otbApplication
#~ import otbApplication
import shapely.geometry
from shapely.geometry import box, Polygon, MultiPolygon, GeometryCollection, mapping, shape
from shapely.ops import transform
import rasterio
from rasterio.features import shapes, rasterize
import numpy as np
from scipy import ndimage
from typing import Sequence, List, Dict, Union
import gdal
import os
......@@ -264,44 +265,79 @@ def create_cloud_mask_v2(
logger.info("Creating cloud-mask_v2")
out_temp_path = Path(Config().get("temp_path"))
out_bin = str(out_temp_path / (out_path.stem + "tmp_bin.tif"))
out_erode = str(out_temp_path / (out_path.stem + "tmp_erode.tif"))
out_dilate = str(out_temp_path / (out_path.stem + "tmp_dilate.tif"))
# Seuillage du masque nuage Sen2Cor
#~ out_bin = str(out_temp_path / (out_path.stem + "_tmp_bin.tif"))
#~ out_erode = str(out_temp_path / (out_path.stem + "_tmp_erode.tif"))
out_dilate = str(out_temp_path / (out_path.stem + "_tmp_dilate.tif"))
logger.info('Loading cloud_prb...')
CLD_seuil = 25
BandMath = otbApplication.Registry.CreateApplication("BandMath")
BandMath.SetParameterStringList("il", [str(cloud_mask)])
BandMath.SetParameterString("out", out_bin)
BandMath.SetParameterOutputImagePixelType("out", 0)
BandMath.SetParameterString("exp", "im1b1 > " + str(CLD_seuil))
BandMath.ExecuteAndWriteOutput()
# Erosion
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
BinaryMorphologicalOperation.SetParameterString("in", out_bin)
BinaryMorphologicalOperation.SetParameterString("out", out_erode)
BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
BinaryMorphologicalOperation.SetParameterInt("channel", 1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", 1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", 1)
BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter", "erode")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
# Dilatation
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
BinaryMorphologicalOperation.SetParameterString("in", out_erode)
BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
BinaryMorphologicalOperation.SetParameterInt("channel", 1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation)
BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
with rasterio.open(str(cloud_mask)) as cld_src:
cld = cld_src.read(1).astype(np.int16)
cld_profile = cld_src.profile
cld_bin = np.where(cld > CLD_seuil, 1, 0).astype(np.uint8)
cld_erode = ndimage.binary_erosion(cld_bin).astype(cld_bin.dtype)
kernel = np.array([[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0]])
cld_dilate = ndimage.binary_dilation(cld_erode, kernel).astype(cld_erode.dtype)
logger.info("saving TIFF file...")
cld_profile.update(driver="Gtiff",
compress="NONE",
tiled=False,
dtype=np.int8,
transform=cld_src.transform,
count=1)
cld_profile.pop('tiled', None)
cld_profile.pop('nodata', None)
with rasterio.Env(GDAL_CACHEMAX=512) as env:
with rasterio.open(out_dilate, "w", **cld_profile) as dst:
dst.write(cld_dilate.astype(np.int8), 1)
#~ # Seuillage du masque nuage Sen2Cor
#~ CLD_seuil = 25
#~ BandMath = otbApplication.Registry.CreateApplication("BandMath")
#~ BandMath.SetParameterStringList("il", [str(cloud_mask)])
#~ BandMath.SetParameterString("out", out_bin)
#~ BandMath.SetParameterOutputImagePixelType("out", 0)
#~ BandMath.SetParameterString("exp", "im1b1 > " + str(CLD_seuil))
#~ BandMath.ExecuteAndWriteOutput()
#~ # Erosion
#~ BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
#~ BinaryMorphologicalOperation.SetParameterString("in", out_bin)
#~ BinaryMorphologicalOperation.SetParameterString("out", out_erode)
#~ BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
#~ BinaryMorphologicalOperation.SetParameterInt("channel", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", 1)
#~ BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", 1)
#~ BinaryMorphologicalOperation.SetParameterString("filter", "erode")
#~ BinaryMorphologicalOperation.ExecuteAndWriteOutput()
#~ # Dilatation
#~ BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
#~ BinaryMorphologicalOperation.SetParameterString("in", out_erode)
#~ BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
#~ BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
#~ BinaryMorphologicalOperation.SetParameterInt("channel", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation)
#~ BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
#~ BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
#~ BinaryMorphologicalOperation.ExecuteAndWriteOutput()
# Save to JP2000
logger.info("converting to JP2000 file...")
src_ds = gdal.Open(out_dilate)
driver = gdal.GetDriverByName("JP2OpenJPEG")
dst_ds = driver.CreateCopy(str(out_path), src_ds,
......@@ -309,8 +345,8 @@ def create_cloud_mask_v2(
dst_ds = None
src_ds = None
os.remove(out_bin)
os.remove(out_erode)
#~ os.remove(out_bin)
#~ os.remove(out_erode)
os.remove(out_dilate)
......@@ -332,72 +368,114 @@ def create_cloud_mask_b11(
logger.info("Masking cloud-mask_v2 with B11")
out_temp_path = Path(Config().get("temp_path"))
out_bin = str(out_temp_path / (out_path.stem + "_tmp_bin.tif"))
out_dilate = str(out_temp_path / (out_path.stem + "_tmp_dilate.tif"))
#~ out_bin = str(out_temp_path / (out_path.stem + "_tmp_bin.tif"))
#~ out_dilate = str(out_temp_path / (out_path.stem + "_tmp_dilate.tif"))
out_mask = str(out_temp_path / (out_path.stem + "_tmp_mask.tif"))
# Seuillage du masque B1
logger.info('Loading B11...')
b11_seuil = 1500
BandMath = otbApplication.Registry.CreateApplication("BandMath")
BandMath.SetParameterStringList("il", [str(b11_path)])
BandMath.SetParameterString("out", out_bin)
BandMath.SetParameterOutputImagePixelType("out", 0)
BandMath.SetParameterString("exp", "im1b1 < " + str(b11_seuil))
BandMath.ExecuteAndWriteOutput()
#~ # Erosion
with rasterio.open(str(b11_path)) as b11_src:
b11 = b11_src.read(1).astype(np.int16)
b11_bin = np.where(b11 < b11_seuil, 1, 0).astype(np.uint8)
kernel = np.array([[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0]])
b11_dilate = ndimage.binary_dilation(b11_bin, kernel).astype(b11_bin.dtype)
with rasterio.open(str(cloud_mask)) as cld_src:
cld = cld_src.read(1).astype(np.int16)
cld_profile = cld_src.profile
cld_mskd = ((cld == 1) * (b11_dilate == 0)).astype(np.uint8)
kernel = np.array([[0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0]])
cld_mskd_dilate = ndimage.binary_dilation(cld_mskd, kernel).astype(cld_mskd.dtype)
logger.info("saving TIFF file...")
cld_profile.update(driver="Gtiff",
compress="NONE",
tiled=False,
dtype=np.int8,
transform=cld_src.transform,
count=1)
cld_profile.pop('tiled', None)
cld_profile.pop('nodata', None)
with rasterio.Env(GDAL_CACHEMAX=512) as env:
with rasterio.open(out_mask, "w", **cld_profile) as dst:
dst.write(cld_mskd_dilate.astype(np.int8), 1)
#~ # Seuillage du masque B11
#~ b11_seuil = 1500
#~ BandMath = otbApplication.Registry.CreateApplication("BandMath")
#~ BandMath.SetParameterStringList("il", [str(b11_path)])
#~ BandMath.SetParameterString("out", out_bin)
#~ BandMath.SetParameterOutputImagePixelType("out", 0)
#~ BandMath.SetParameterString("exp", "im1b1 < " + str(b11_seuil))
#~ BandMath.ExecuteAndWriteOutput()
#~ # Dilatation
#~ BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
#~ BinaryMorphologicalOperation.SetParameterString("in", out_bin)
#~ BinaryMorphologicalOperation.SetParameterString("out", out_erode)
#~ BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
#~ BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
#~ BinaryMorphologicalOperation.SetParameterInt("channel", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", 1)
#~ BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", 1)
#~ BinaryMorphologicalOperation.SetParameterString("filter", "erode")
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation)
#~ BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
#~ BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
#~ BinaryMorphologicalOperation.ExecuteAndWriteOutput()
# Dilatation
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
BinaryMorphologicalOperation.SetParameterString("in", out_bin)
BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
BinaryMorphologicalOperation.SetParameterInt("channel", 1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation)
BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
# Masquage du masque
BandMath = otbApplication.Registry.CreateApplication("BandMath")
BandMath.SetParameterStringList("il", [str(cloud_mask),str(out_dilate)])
BandMath.SetParameterString("out", out_mask)
BandMath.SetParameterOutputImagePixelType("out", 0)
BandMath.SetParameterString("exp", "(im1b1 == 1) * (im2b1 == 0)")
BandMath.ExecuteAndWriteOutput()
# Dilatation
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
BinaryMorphologicalOperation.SetParameterString("in", out_mask)
BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
BinaryMorphologicalOperation.SetParameterInt("channel", 1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation+1)
BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation+1)
BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
#~ BandMath = otbApplication.Registry.CreateApplication("BandMath")
#~ BandMath.SetParameterStringList("il", [str(cloud_mask),str(out_dilate)])
#~ BandMath.SetParameterString("out", out_mask)
#~ BandMath.SetParameterOutputImagePixelType("out", 0)
#~ BandMath.SetParameterString("exp", "(im1b1 == 1) * (im2b1 == 0)")
#~ BandMath.ExecuteAndWriteOutput()
#~ # Dilatation
#~ BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
#~ BinaryMorphologicalOperation.SetParameterString("in", out_mask)
#~ BinaryMorphologicalOperation.SetParameterString("out", out_dilate)
#~ BinaryMorphologicalOperation.SetParameterOutputImagePixelType("out", 0)
#~ BinaryMorphologicalOperation.SetParameterInt("channel", 1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.xradius", dilatation+1)
#~ BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", dilatation+1)
#~ BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
#~ BinaryMorphologicalOperation.SetParameterString("filter", "dilate")
#~ BinaryMorphologicalOperation.ExecuteAndWriteOutput()
# Save to JP2000
src_ds = gdal.Open(out_dilate)
logger.info("converting to JP2000 file...")
src_ds = gdal.Open(out_mask)
driver = gdal.GetDriverByName("JP2OpenJPEG")
dst_ds = driver.CreateCopy(str(out_path), src_ds,
options=["CODEC=JP2", "QUALITY=100", "REVERSIBLE=YES", "YCBCR420=NO"])
dst_ds = None
src_ds = None
os.remove(out_bin)
os.remove(out_dilate)
#~ os.remove(out_bin)
#~ os.remove(out_dilate)
os.remove(out_mask)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment