diff --git a/TimeSeries/landsat_planetary.py b/TimeSeries/landsat_planetary.py
index aa2db6a897a9ce51c8313ce629efd461e37b1c44..881d6cf34a33d4b7ef22e1c99db266e5298b4ab8 100644
--- a/TimeSeries/landsat_planetary.py
+++ b/TimeSeries/landsat_planetary.py
@@ -1,3 +1,4 @@
+from TimeSeries.s2theia import *
from Common.geotools import get_query_bbox
from pystac_client import Client
import planetary_computer as PC
@@ -52,4 +53,133 @@ def fetch(shp, dt, out_fld, auth=None, band_list=None):
with rasterio.open(ofn, "w", **out_meta) as dest:
dest.write(out_img)
prg.update()
- prg.close()
\ No newline at end of file
+ prg.close()
+
+
+class Landsat89PlanetaryTilePipeline(S2TheiaTilePipeline):
+ # --- BEGIN SENSOR PROTOTYPE ---
+
+ NAME = 'LANDSAT89-PLANETARY'
+ REF_TYPE = otb.ImagePixelType_uint16
+ MSK_TYPE = otb.ImagePixelType_uint8
+ PTRN_dir = 'LC0*_L2SP*'
+ PTRN_ref = '_B'
+ B2_name = 'B1'
+ PTRN_30m = ['*_SR_B1.TIF', '*_SR_B2.TIF', '*_SR_B3.TIF',
+ '*_SR_B4.TIF', '*_SR_B5.TIF', '*_SR_B6.TIF',
+ '*_SR_B7.TIF']
+ PTRN_msk = ['*_QA_PIXEL.TIF']
+ PTRN_ful = PTRN_30m
+ MASK_exp = 'im1b1!=21824 && im1b1!=21888 && im1b1!=30048'
+ FEAT_exp = {
+ 'B1': 'im1b1',
+ 'B2': 'im1b2',
+ 'B3': 'im1b3',
+ 'B4': 'im1b4',
+ 'B5': 'im1b5',
+ 'B6': 'im1b6',
+ 'B7': 'im1b7',
+ 'NDVI': '(im1b5-im1b4)/(im1b5+im1b4+1e-6)',
+ 'NDWI': '(im1b3-im1b5)/(im1b3+im1b5+1e-6)',
+ 'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 8)]) + ')',
+ 'MNDWI': '(im1b3-im1b7)/(im1b3+im1b7+1e-6)',
+ 'SWNDVI': '(im1b7-im1b5)/(im1b7+im1b5+1e-6)',
+ }
+ NDT = 0
+
+ @classmethod
+ def _check(cls,x):
+ return os.path.isdir(x) and os.path.basename(x).startswith(cls.PTRN_dir.split('*')[0]) and 'L2SP' in os.path.basename(x)
+
+ @classmethod
+ def _img_id(cls,x):
+ if cls._check(x):
+ return os.path.basename(x)
+ else:
+ return None
+
+ @classmethod
+ def _img_date(cls,x):
+ if cls._check(x):
+ return os.path.basename(x).split('_')[3]
+ else:
+ return None
+
+ @classmethod
+ def _tile_id(cls,x):
+ if cls._check(x):
+ return os.path.basename(x).split('_')[2]
+ else:
+ return None
+
+ def _process_mask(self, msks):
+ msk_pipe = otb.Registry.CreateApplication('BandMath')
+ msk_pipe.AddImageToParameterInputImageList('il', msks[0].GetParameterOutputImage('out'))
+ msk_pipe.SetParameterString('exp', self.MASK_exp)
+ msk_pipe.Execute()
+ return msk_pipe
+
+ # ---- END SENSOR PROTOTYPE ----
+
+ def preprocess(self):
+ proc_idx = self.out_idx.copy()
+ self.out_idx = []
+ num_of_band_inputs = len(self.PTRN_ful)
+ num_of_mask_inputs = len(self.PTRN_msk)
+ for t in proc_idx[::num_of_band_inputs+num_of_mask_inputs]:
+ idx_to_stack = []
+ cct = otb.Registry.CreateApplication('ConcatenateImages')
+ [cct.AddImageToParameterInputImageList('il', self.pipe[j].GetParameterOutputImage('out')) for j in range(t, t+num_of_band_inputs)]
+ cct.Execute()
+ fn = self.files[t].replace(self.B2_name, 'STACK')
+ ty = self.REF_TYPE
+ self.append(cct, fn, ty, 'out', is_output=True)
+
+ msks = self.pipe[t+num_of_band_inputs:t+num_of_band_inputs+num_of_mask_inputs]
+ msk_pipe = self._process_mask(msks)
+ fn = self.files[t].replace(self.B2_name, 'BINARY_MASK')
+ ty = self.MSK_TYPE
+ self.append(msk_pipe, fn, ty, 'out', is_output=True)
+
+class Landsat7PlanetaryTilePipeline(Landsat89PlanetaryTilePipeline):
+ NAME = 'LANDSAT7-PLANETARY'
+ PTRN_dir = 'LE07_L2SP*'
+ PTRN_30m = ['*_SR_B1.TIF', '*_SR_B2.TIF', '*_SR_B3.TIF',
+ '*_SR_B4.TIF', '*_SR_B5.TIF', '*_SR_B7.TIF']
+ PTRN_ful = PTRN_30m
+ MASK_exp = 'im1b1!=5440 && im1b1!=5504 && im1b1!=13664'
+ FEAT_exp = {
+ 'B1': 'im1b1',
+ 'B2': 'im1b2',
+ 'B3': 'im1b3',
+ 'B4': 'im1b4',
+ 'B5': 'im1b5',
+ 'B7': 'im1b6',
+ 'NDVI': '(im1b4-im1b3)/(im1b4+im1b3+1e-6)',
+ 'NDWI': '(im1b2-im1b4)/(im1b2+im1b4+1e-6)',
+ 'BRI': 'sqrt(' + '+'.join(['im1b%d*im1b%d' % (i, i) for i in range(1, 7)]) + ')',
+ 'MNDWI': '(im1b2-im1b6)/(im1b2+im1b6+1e-6)',
+ 'SWNDVI': '(im1b6-im1b4)/(im1b6+im1b4+1e-6)',
+ }
+
+class Landsat45PlanetaryTilePipeline(Landsat7PlanetaryTilePipeline):
+ NAME = 'LANDSAT45-PLANETARY'
+ PTRN_dir = 'LT0*_L2SP*'
+
+class Landsat89PlanetaryPipeline(S2TheiaPipeline):
+
+ S2TilePipeline = Landsat89PlanetaryTilePipeline
+ _check = S2TilePipeline._check
+ _tile_id = S2TilePipeline._tile_id
+
+class Landsat7PlanetaryPipeline(S2TheiaPipeline):
+
+ S2TilePipeline = Landsat7PlanetaryTilePipeline
+ _check = S2TilePipeline._check
+ _tile_id = S2TilePipeline._tile_id
+
+class Landsat45PlanetaryPipeline(S2TheiaPipeline):
+
+ S2TilePipeline = Landsat45PlanetaryTilePipeline
+ _check = S2TilePipeline._check
+ _tile_id = S2TilePipeline._tile_id
diff --git a/TimeSeries/s2theia.py b/TimeSeries/s2theia.py
index b42992cbcda0afb50b79ba6ec1cfbaa4de15a60b..0fec83dd5c6391cce096312653acb8321a8ae788 100644
--- a/TimeSeries/s2theia.py
+++ b/TimeSeries/s2theia.py
@@ -163,9 +163,9 @@ class S2TheiaTilePipeline:
def _check_roi(cls, x, roi, min_surf=0.0, temp_fld='/tmp'):
if cls._check(x):
if os.path.isdir(x):
- bnd = glob.glob(os.path.join(x, cls.PTRN_20m[0]))[0]
+ bnd = glob.glob(os.path.join(x, cls.PTRN_ful[0]))[0]
elif os.path.splitext(x)[-1] == '.zip':
- idf = cls.PTRN_20m[0].replace('*', '')
+ idf = cls.PTRN_ful[0].replace('*', '')
arch = zipfile.ZipFile(x)
fid = [name for name in arch.namelist() if idf in name]
fle = arch.read(fid[0])
@@ -340,7 +340,7 @@ class S2TheiaTilePipeline:
i += 1
to_merge.append(curr)
- T = 10 + len(self.PTRN_msk)
+ T = len(self.PTRN_ful) + len(self.PTRN_msk)
new_dates = []
proc_idx = self.out_idx.copy()
self.out_idx = []
@@ -351,18 +351,18 @@ class S2TheiaTilePipeline:
idx = T*l[0]+k
self.append(self.pipe[idx], self.files[idx], self.types[idx], 'out', is_output=True)
else:
- for k in range(10):
+ for k in range(len(self.PTRN_ful)):
mos = otb.Registry.CreateApplication('Mosaic')
for i in l:
mos.AddImageToParameterInputImageList('il', self.pipe[T*i+k].GetParameterOutputImage(self.out_p[T*i+k]))
mos.SetParameterInt('nodata', self.NDT)
mos.Execute()
self.append(mos, self.files[T*i+k], self.types[T*i+k], 'out', is_output=True)
- for k in range(10,T):
+ for k in range(len(self.PTRN_ful),T):
bm = otb.Registry.CreateApplication('BandMath')
for i in l:
bm.AddImageToParameterInputImageList('il', self.pipe[T*i+k].GetParameterOutputImage(self.out_p[T*i+k]))
- bm.SetParameterString('exp', self.MERG_msk[k-10] + '({})'.format(','.join(['im{}b1'.format(w+1) for w in range(len(l))])))
+ bm.SetParameterString('exp', self.MERG_msk[k-len(self.PTRN_ful)] + '({})'.format(','.join(['im{}b1'.format(w+1) for w in range(len(l))])))
bm.Execute()
self.append(bm, self.files[T * i + k], self.types[T * i + k], 'out', is_output=True)
new_dates.append(self.input_dates[l[-1]])
@@ -385,7 +385,7 @@ class S2TheiaTilePipeline:
return cc
def set_input_epsg(self):
- f, _ = self.get_file(self.image_list[0], self.PTRN_20m[0])
+ f, _ = self.get_file(self.image_list[0], self.PTRN_ful[0])
ds = gdal.Open(f)
self.input_epsg = osr.SpatialReference(wkt=ds.GetProjection()).GetAuthorityCode('PROJCS')
ds = None
@@ -417,10 +417,10 @@ class S2TheiaTilePipeline:
proc_idx = self.out_idx.copy()
self.out_idx = []
num_of_mask_inputs = len(self.PTRN_msk)
- for t in proc_idx[::10+num_of_mask_inputs]:
+ for t in proc_idx[::len(self.PTRN_ful)+num_of_mask_inputs]:
idx_to_stack = []
- for k in range(t,t+10):
- if k % (10 + num_of_mask_inputs) in [3, 4, 5, 7, 8, 9]:
+ for k in range(t,t+len(self.PTRN_ful)):
+ if k % (len(self.PTRN_ful) + num_of_mask_inputs) in [3, 4, 5, 7, 8, 9]:
si = otb.Registry.CreateApplication('Superimpose')
si.SetParameterInputImage('inm', self.pipe[k].GetParameterOutputImage('out'))
si.SetParameterInputImage('inr', self.pipe[t].GetParameterOutputImage('out'))
@@ -442,7 +442,7 @@ class S2TheiaTilePipeline:
ty = self.REF_TYPE
self.append(cct, fn, ty, 'out', is_output=True)
- msks = self.pipe[t+10:t+10+num_of_mask_inputs]
+ msks = self.pipe[t+len(self.PTRN_ful):t+len(self.PTRN_ful)+num_of_mask_inputs]
msk_pipe = self._process_mask(msks)
for app in msk_pipe[:-1]:
self.append(app)
@@ -554,8 +554,8 @@ class S2TheiaTilePipeline:
proc_idx = self.out_idx.copy()
self.out_idx = []
img = [self.pipe[t] for t in range(self.pipe_start+match_band,
- self.pipe_start+(10+len(self.PTRN_msk))*len(self.input_dates),
- 10+len(self.PTRN_msk))]
+ self.pipe_start+(len(self.PTRN_ful)+len(self.PTRN_msk))*len(self.input_dates),
+ len(self.PTRN_ful)+len(self.PTRN_msk))]
msk = [self.pipe[t] for t in proc_idx[1::2]]
RX, RY = img[0].GetImageSpacing(out_param)
ref_idx = None
@@ -662,7 +662,7 @@ class S2TheiaTilePipeline:
gf = otb.Registry.CreateApplication('ImageTimeSeriesGapFilling')
gf.SetParameterInputImage('in', self.pipe[-2].GetParameterOutputImage('out'))
gf.SetParameterInputImage('mask', self.pipe[-1].GetParameterOutputImage('out'))
- gf.SetParameterInt('comp', 10)
+ gf.SetParameterInt('comp', len(self.PTRN_ful))
gf.SetParameterString('it', 'linear')
gf.SetParameterString('id', os.path.join(self.folder, self.tile_id + '_indates.txt'))
gf.SetParameterString('od', output_dates)
@@ -670,7 +670,7 @@ class S2TheiaTilePipeline:
gf.Execute()
gf_out = gf
else:
- gf_fn = self.temp_fld + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILLED_FULL.tif'
+ gf_fn = self.temp_fld + os.sep + '{}_'.format(self.NAME) + self.tile_id + '_GAPFILLED_FULL.tif'
if not os.path.exists(gf_fn):
gf.SetParameterString("out", gf_fn)
gf.ExecuteAndWriteOutput()
@@ -679,14 +679,14 @@ class S2TheiaTilePipeline:
t = 1
for d in od:
- ch_list = ['Channel%d' % i for i in range(t,t+10)]
- t += 10
+ ch_list = ['Channel%d' % i for i in range(t,t+len(self.PTRN_ful))]
+ t += len(self.PTRN_ful)
er = otb.Registry.CreateApplication('ExtractROI')
er.SetParameterInputImage('in', gf_out.GetParameterOutputImage('out'))
er.UpdateParameters()
er.SetParameterStringList('cl', ch_list)
er.Execute()
- dn = 'SENTINEL2_' + self.tile_id + '_GAPFILL_' + d
+ dn = '{}_'.format(self.NAME) + self.tile_id + '_GAPFILL_' + d
fn = self.tile_id + os.sep + dn + os.sep + dn + '_STACK.tif'
ty = self.REF_TYPE
self.append(er, fn, ty, 'out', is_output=True)
@@ -726,7 +726,7 @@ class S2TheiaTilePipeline:
i += 1
bm.SetParameterString('exp', '{' + ';'.join(expr) + '}')
bm.Execute()
- fn = self.tile_id + os.sep + 'SENTINEL2_' + self.tile_id + '_GAPFILL_' + feat + '.tif'
+ fn = self.tile_id + os.sep + '{}_'.format(self.NAME) + self.tile_id + '_GAPFILL_' + feat + '.tif'
ty = otb.ImagePixelType_float
self.append(bm, fn, ty, 'out', is_output=True)
@@ -892,7 +892,7 @@ class S2TheiaPipeline:
out_mos = []
vrtopt = gdal.BuildVRTOptions()
for i in range(n_dates):
- fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + self.tiles[0].output_dates[i] + '_' + stack_name + '.vrt'
+ fn = out_fld + os.sep + '{}_MOSAIC_GAPFILL_'.format(self.tiles[0].NAME) + self.tiles[0].output_dates[i] + '_' + stack_name + '.vrt'
to_mosaic = [x[i] for x in out]
gdal.BuildVRT(fn, to_mosaic, options=vrtopt)
out_mos.append(fn)
@@ -916,7 +916,7 @@ class S2TheiaPipeline:
out_mos = []
vrtopt = gdal.BuildVRTOptions()
for i in range(len(feat_list)):
- fn = out_fld + os.sep + 'SENTINEL2_MOSAIC_GAPFILL_' + feat_list[i] + '.vrt'
+ fn = out_fld + os.sep + '{}_MOSAIC_GAPFILL_'.format(self.tiles[0].NAME) + feat_list[i] + '.vrt'
to_mosaic = [x[i] for x in out]
gdal.BuildVRT(fn, to_mosaic, options=vrtopt)
out_mos.append(fn)
diff --git a/Workflows/basic.py b/Workflows/basic.py
index 81fd68242d69630a1669b254e406b8b85cb66dda..6d439b86e5fb4db50f5c107735d8daacc6f221b6 100644
--- a/Workflows/basic.py
+++ b/Workflows/basic.py
@@ -51,6 +51,16 @@ def process_timeseries(oroot, d, ts_lst_pkl):
os.makedirs(ots, exist_ok=True)
ts_lst.append(preprocess_planet(ts['path'],
ots))
+ elif ts['type'].startswith('landsat'):
+ print('[MORINGA-INFO] : Preprocessing {}'.format(ts['type']))
+ ots = os.path.join(oroot, 'timeseries/' + ts['type'])
+ os.makedirs(ots, exist_ok=True)
+ ts_lst.append(preprocess_landsat(ts['path'],
+ ots,
+ roi=d['roi'],
+ output_dates_file=ts['output_dates_file'],
+ which=ts['type']))
+
else:
raise ValueError('TimeSeries type not yet supported.')
with open(ts_lst_pkl, 'wb') as ts_save:
@@ -84,7 +94,9 @@ def train_valid_workflow(seg, ts_lst_pkl, d, m_file):
reference_data=d['ref_db']['path'],
ref_class_field=d['ref_db']['fields'])
- obc.gen_k_folds(5, class_field=d['ref_db']['fields'][-1],augment=d['training']['augment_if_missing'])
+ obc.gen_k_folds(5, class_field=d['ref_db']['fields'][-1],
+ augment=d['training']['augment_if_missing'],
+ min_samples_per_class=d['training']['augment_if_missing'])
if 'export_training_base' in d['training'].keys() and d['training']['export_training_base'] is True:
obc.save_training_base('{}/_side/training_base.pkl'.format(os.path.join(d['output_path'], d['chain_name'])))
diff --git a/Workflows/operations.py b/Workflows/operations.py
index 3619f5f40da921dac5c82706c47afafd491e35f5..a80026c680a3d2d9620876e50f9d95c59332055e 100644
--- a/Workflows/operations.py
+++ b/Workflows/operations.py
@@ -55,6 +55,25 @@ def preprocess_s2(in_fld, out_fld, output_dates_file=None, roi=None,
return s2.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt', align=align,
align_to=align_to, align_to_band=align_to_band, align_using_band=align_using_band)
+def preprocess_landsat(in_fld, out_fld, output_dates_file=None, roi=None, which='landsat89'):
+ LandsatProcessor = None
+ if which == 'landsat89':
+ LandsatProcessor = landsat_planetary.Landsat89PlanetaryPipeline
+ elif which == 'landsat7':
+ LandsatProcessor = landsat_planetary.Landsat7PlanetaryPipeline
+ elif which == 'landsat45':
+ LandsatProcessor = landsat_planetary.Landsat45PlanetaryPipeline
+ else:
+ raise ValueError("Unknown Landsat mission.")
+
+ ls = LandsatProcessor(in_fld, temp_fld=out_fld, roi=roi)
+ if output_dates_file is not None:
+ ls.set_output_dates_by_file(output_dates_file)
+ else:
+ raise ValueError("Please provide path to a text file containing output dates.")
+
+ return ls.extract_feature_set(out_fld, store_gapfill=True, mosaicking='vrt')
+
def preprocess_planet(in_fld, out_fld):
pl = planet_mosaics.PlanetMosaicPipeline(in_fld)
pl.compute_features()