DTACSNet: onboard cloud detection and atmospheric correction with end-to-end deep learning emulators
This repo contains an open implementation to run inference with DTACSNet models for atmospheric correction and also for cloud detection. These two models are independent (you could run one or the other or both). The trained models provided here are customized to the band configuration that will be available in Phi-Sat-II (excluding PAN band). Trained models are released under a Creative Commons non-commercial licence
.
Phi-Sat-II simulated image run
In this notebook we run the DTACSNet model over a simulated Phi-Sat-II image. To run this tutorial we downloaded this image from here and save it in the tutorials/examples folder.
# download the image
import os
if not os.path.exists("examples/phisat2-l1c.tiff"):
!wget -P examples/ https://cloud.sinergise.com/s/RSgHGwWWCxYp5yq/download/phisat2-l1c.tiff
import sys
sys.path.append("..")
Load data
Load the simulated Phi-Sat-II image. We multiply the values by 10,000 to follow the same convention as in Sentinel-2.
import rasterio
import rasterio.windows
import os
phisat2_bands = ["B2","B3","B4","B5","B6","B7","B8"]
PHISAT2_PAN_BANDS = ["B2", "B3", "B4", "PAN", "B8", "B5", "B6", "B7"]
folder_examples = "examples"
with rasterio.open(os.path.join(folder_examples,"phisat2-l1c.tiff")) as rst:
indexes_read = [PHISAT2_PAN_BANDS.index(b) + 1 for b in phisat2_bands]
data = rst.read(indexes_read).astype("float32") * 10_000 # Convert to the same units as S2 TOA images
transform_phisat2 = rst.transform
data.shape
# Data are TOA reflectance values. We multiply them by 10_000 to use the same convention as in Sentinel-2
data
Atmospheric correction model
Load model
from dtacs.model_wrapper import ACModel
import torch
model_atmospheric_correction = ACModel(model_name="CNN_corrector_phisat2")
model_atmospheric_correction.load_weights()
Run inference
%%time
ac_output = model_atmospheric_correction.predict(data)
ac_output
Plot
import matplotlib.pyplot as plt
from dtacs import plot
import rasterio.plot as rstplt
fig, ax = plt.subplots(2,2,figsize=(12,12),tight_layout=True)
rstplt.show(data[2::-1,...]/4_000,ax=ax[0,0])
ax[0,0].set_title("RGB PhiSat-II L1C")
rstplt.show(ac_output[2::-1,...]/4_000,ax=ax[0,1])
ax[0,1].set_title("RGB PhiSat-II L2A (DTACSNet corrected)")
nirredgreen = [-1, 2, 1]
rstplt.show(data[nirredgreen,...]/10_000,ax=ax[1,0])
ax[1,0].set_title("NIRRG PhiSat-II L1C")
rstplt.show(ac_output[nirredgreen,...]/10_000,ax=ax[1,1])
ax[1,1].set_title("NIRRG PhiSat-II L2A (DTACSNet corrected)")
Compare L2A Phi-Sat-II image with L2A Sentinel-2 image
The simulated Phi-Sat-II image that we have shown below has been created from the Sentinel-2 tile: S2A_MSIL1C_20191111T132241_N0208_R038_T23LLK_20191111T145626. We have copied a 500x500 subset of this tile in the tutorials/examples folder. We have also downloaded the corresponding L2A Sentinel-2 file to compare the atmospheric correction output of the two products.
# Load Sentinel-2 L2A
with rasterio.open(os.path.join(folder_examples,"S2L2A.tif")) as rst:
indexes_read = [rst.descriptions.index(b) + 1 for b in phisat2_bands]
data_l2a = rst.read(indexes_read)
# Load Sentinel-2 L1C
with rasterio.open(os.path.join(folder_examples,"S2L1C.tif")) as rst:
indexes_read = [rst.descriptions.index(b) + 1 for b in phisat2_bands]
data_l1c = rst.read(indexes_read)
bounds = rst.bounds
window_in = rasterio.windows.from_bounds(*bounds, transform=transform_phisat2)
slice_read = (slice(None),) + window_in.toslices()
fig, ax = plt.subplots(2,4,figsize=(24,12),tight_layout=True)
rstplt.show(data[2::-1,...][slice_read]/4_000,ax=ax[0,0])
ax[0,0].set_title("RGB PhiSat-II L1C")
rstplt.show(ac_output[2::-1,...][slice_read]/4_000,ax=ax[0,1])
ax[0,1].set_title("RGB PhiSat-II L2A (DTACSNet corrected)")
rstplt.show(data_l1c[2::-1,...]/4_000,ax=ax[0,2])
ax[0,2].set_title("RGB Sentinel-2 L1C")
rstplt.show(data_l2a[2::-1,...]/4_000,ax=ax[0,3])
ax[0,3].set_title("RGB Sentinel-2 L2A")
nirredgreen = [-1, 2, 1]
rstplt.show(data[nirredgreen,...][slice_read]/10_000,ax=ax[1,0])
ax[1,0].set_title("NIRRG PhiSat-II L1C")
rstplt.show(ac_output[nirredgreen,...][slice_read]/10_000,ax=ax[1,1])
ax[1,1].set_title("NIRRG PhiSat-II L2A (DTACSNet corrected)")
rstplt.show(data_l1c[nirredgreen,...]/10_000,ax=ax[1,2])
ax[1,2].set_title("NIRRG PhiSat-II L1C")
rstplt.show(data_l2a[nirredgreen,...]/10_000,ax=ax[1,3])
ax[1,3].set_title("NIRRG Sentinel-2 L2A")
