Temporal Resampling and gapfilling
The package
otbApplicationis not installed in SISE environment. You must activate SISE-3.7 environnment by launching the following command in Anaconda Prompt :conda activate sise-3.7
[16]:
import glob, os, subprocess
import pandas as pd
import otbApplication as otb
from pathlib import Path
import datetime
import rasterio
import numpy as np
print('All libraries successfully imported!')
All libraries successfully imported!
Set directory
[17]:
computer_path = '/export/miro/ndeffense/LBRAT2104/'
grp_nb = '99'
# Directory for all work files
work_path = f'{computer_path}STUDENTS/GROUP_{grp_nb}/TP/'
clipped_path = f'{work_path}2_L2A_CLIPPED/'
gapfilled_path = f'{work_path}4_L2A_GAPFILLED/'
Path(gapfilled_path).mkdir(parents=True, exist_ok=True)
Parameters
[18]:
band_list = ['B02','B03']
format = '%Y%m%d' # format : YYYYMMDD
begin_date = 7
# Distance between interpolated dates
frequency = '1m'
# Interpolation type
interp_type = 'linear' # 'spline'
flag_otbAPI = False
Prepare data
Get dataframe with all inputs
[19]:
dict_list = []
for im in glob.glob(f'{clipped_path}*SCL*10m*.tif'):
date_str = os.path.basename(im)[begin_date:begin_date+8]
date_im = datetime.datetime.strptime(date_str, format).date()
dict_list.append({'date': date_im,
'date_str': date_str,
'SCL': im})
df = pd.DataFrame.from_dict(dict_list).sort_values('date')
for band in band_list:
df[band] = "" # Create empty column in dataframe
for im in glob.glob(f'{clipped_path}*{band}*.tif'):
date_str = os.path.basename(im)[begin_date:begin_date+8]
date_im = datetime.datetime.strptime(date_str, format).date()
df.loc[df.date == date_im, band] = im
display(df)
| date | date_str | SCL | B02 | B03 | |
|---|---|---|---|---|---|
| 0 | 2020-01-16 | 20200116 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 1 | 2020-02-12 | 20200212 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 2 | 2020-03-16 | 20200316 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 3 | 2020-04-17 | 20200417 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 4 | 2020-05-20 | 20200520 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 5 | 2020-06-21 | 20200621 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 6 | 2020-07-19 | 20200719 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 7 | 2020-08-13 | 20200813 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 8 | 2020-09-14 | 20200914 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 9 | 2020-10-19 | 20201019 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 10 | 2020-11-18 | 20201118 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
| 11 | 2020-12-18 | 20201218 | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... | /export/miro/ndeffense/LBRAT2104/STUDENTS/GROU... |
Modify and concatenate masks
[20]:
mask_concat_file = f'{gapfilled_path}SCL_timeserie.tif'
list_src_arr = []
for i, row in df.iterrows():
im_file = row['SCL']
src = rasterio.open(im_file, "r")
out_meta = src.meta
out_meta.update(count=len(df))
# Read file as numpy array
SCL = src.read(1)
src.close()
SCL[SCL == 0] = 1 # No data
SCL[SCL == 1] = 1 # Saturated or defective
SCL[SCL == 2] = 1 # Dark area pixels
SCL[SCL == 3] = 1 # Cloud shadows
SCL[SCL == 4] = 0 # Vegetation
SCL[SCL == 5] = 0 # Not vegetated
SCL[SCL == 6] = 0 # Water
SCL[SCL == 7] = 0 # Unclassified
SCL[SCL == 8] = 1 # Cloud medium probability
SCL[SCL == 9] = 1 # Cloud high probability
SCL[SCL == 10] = 1 # Thin cirrus
SCL[SCL == 11] = 1 # Snow
list_src_arr.append(SCL)
mask_arr_stack = np.dstack(list_src_arr)#.astype(np.float32)
print(mask_arr_stack.shape)
print(f'There are {mask_arr_stack.shape[2]} masks')
with rasterio.open(mask_concat_file, "w", **out_meta) as dest:
for band_nr, src in enumerate(list_src_arr, start=1):
dest.write(src, band_nr)
(570, 986, 12)
There are 12 masks
Get input dates
[21]:
input_dates_txt = f'{gapfilled_path}input_dates.txt'
input_dates_list = df['date_str'].to_list()
with open(input_dates_txt, 'w') as f:
for i in range(0,len(input_dates_list)):
f.write(input_dates_list[i] + '\n')
Get output dates
[22]:
output_dates_txt = f'{gapfilled_path}output_dates.txt'
start_date = df['date'].iloc[0]
end_date = df['date'].iloc[-1]
output_dates_list = pd.date_range(start_date,end_date,freq=frequency).strftime('%Y%m%d').to_list()
with open(output_dates_txt, 'w') as f:
for i in range(0,len(output_dates_list)):
print(output_dates_list[i])
f.write(output_dates_list[i] + '\n')
20200131
20200229
20200331
20200430
20200531
20200630
20200731
20200831
20200930
20201031
20201130
Concatenate reflectance timeseries
[27]:
for band in band_list:
im_concat_file = f'{gapfilled_path}{band}_timeserie.tif'
im_path_list = df[band].tolist()
if flag_otbAPI:
app = otb.Registry.CreateApplication("ConcatenateImages")
app.SetParameterStringList("il", im_path_list)
app.SetParameterString("out", im_concat_file)
app.ExecuteAndWriteOutput()
else:
cmd = 'otbcli_ConcatenateImages'
cmd += f' -il {" ".join(im_path_list)}'
cmd += f' -out {im_concat_file}'
print(cmd)
subprocess.call(cmd, shell=True)
print('-------------------------------------------')
otbcli_ConcatenateImages -il /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200116T105309_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200212T104049_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200316T104709_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200417T104021_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200520T105031_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200621T103629_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200719T105031_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200813T104629_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200914T104031_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201019T103959_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201118T104329_B02_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201218T104349_B02_10m_ROI.tif -out /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B02_timeserie.tif
-------------------------------------------
otbcli_ConcatenateImages -il /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200116T105309_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200212T104049_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200316T104709_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200417T104021_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200520T105031_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200621T103629_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200719T105031_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200813T104629_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20200914T104031_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201019T103959_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201118T104329_B03_10m_ROI.tif /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/2_L2A_CLIPPED/T31UFS_20201218T104349_B03_10m_ROI.tif -out /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B03_timeserie.tif
-------------------------------------------
Temporal Resampling and Gapfilling
[26]:
for band in band_list:
im_concat_file = f'{gapfilled_path}{band}_timeserie.tif'
gapfilled_file = f'{gapfilled_path}{band}_timeserie_gapfilled_{frequency}.tif'
if flag_otbAPI:
app = otb.Registry.CreateApplication("ImageTimeSeriesGapFilling")
app.SetParameterString("in", im_concat_file)
app.SetParameterString("mask", mask_concat_file)
app.SetParameterString("out", gapfilled_file)
app.SetParameterInt("comp", 1)
app.SetParameterString("it", interp_type)
app.SetParameterString("id", input_dates_txt)
app.SetParameterString("od", output_dates_txt)
app.ExecuteAndWriteOutput()
else:
cmd = 'otbcli_ImageTimeSeriesGapFilling'
cmd += f' -in {im_concat_file}'
cmd += f' -mask {mask_concat_file}'
cmd += f' -out {gapfilled_file}'
cmd += f' -comp 1'
cmd += f' -it {interp_type}'
cmd += f' -id {input_dates_txt}'
cmd += f' -od {output_dates_txt}'
print(cmd)
subprocess.call(cmd, shell=True)
print('-------------------------------------------')
otbcli_ImageTimeSeriesGapFilling -in /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B02_timeserie.tif -mask /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/SCL_timeserie.tif -out /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B02_timeserie_gapfilled_1m.tif -comp 1 -it linear -id /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/input_dates.txt -od /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/output_dates.txt
-------------------------------------------
otbcli_ImageTimeSeriesGapFilling -in /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B03_timeserie.tif -mask /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/SCL_timeserie.tif -out /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/B03_timeserie_gapfilled_1m.tif -comp 1 -it linear -id /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/input_dates.txt -od /export/miro/ndeffense/LBRAT2104/STUDENTS/GROUP_99/TP/4_L2A_GAPFILLED/output_dates.txt
-------------------------------------------
Display gapfilled timeserie
TO DO