Commit 0e62de83 authored by Gaetano Raffaele's avatar Gaetano Raffaele
Browse files

WIP: complete refactoring of stratification application.

parent 0465cb36
import gdal, ogr
import datetime
import otbApplication as otb
import os
def find_next_date(datelist, date):
for i,l in enumerate(datelist):
if l >= date:
return i
return None
def find_first_month_day(dates, month=1, day=1):
start_year = dates[0].year
test_date = datetime.datetime(start_year,month,day)
if dates[0] > test_date:
start_year += 1
test_date = datetime.datetime(start_year, month, day)
return find_next_date(dates, test_date)
def get_period_intervals(date_file, md=[1, 1], duration=365):
dates = []
with open(date_file) as f:
for l in f.read().splitlines():
try:
dates.append(datetime.datetime.strptime(l, '%Y%m%d'))
except:
continue
delta = dates[-1] - dates[-2]
dates.append(dates[-1] + delta)
periods = []
s, e = 0, 0
S = 0
while s is not None and e is not None:
s = find_first_month_day(dates, md[0], md[1])
if s is not None:
e = find_next_date(dates, dates[s] + datetime.timedelta(days=duration))
if e is not None:
periods.append((S + s, S + e))
S += s + 1
dates = dates[s + 1:]
return periods
def get_mean_expression(periods):
l = periods[0][1] - periods[0][0] + 1
expr = []
for i in range(1, l + 1):
s = [p[0] + i for p in periods]
expr.append('(' + '+'.join(['im1b%d' % x for x in s]) + ')/' + str(len(s)))
return '{' + ';'.join(expr) + '}'
def temporal_smoothing(file_list, output_file):
nd = len(file_list)
ds = gdal.Open(file_list[0])
ndv = ds.GetRasterBand(1).GetNoDataValue()
ds = None
dates = [x.split('_')[-2][3:] for x in file_list]
sdates = [datetime.datetime.strptime(d, '%Y%j').date().strftime('%Y%m%d') for d in dates]
datefile = os.path.splitext(output_file)[0] + '_dates.txt'
with open(datefile, 'w') as df:
df.write('\n'.join(sdates) + '\n')
sg = otb.Registry.CreateApplication('TemporalSmoothing')
sg.SetParameterStringList('il', file_list)
sg.SetParameterString('dates', datefile)
sg.Execute()
msks = otb.Registry.CreateApplication('BandMathX')
msks.AddImageToParameterInputImageList('il', sg.GetParameterOutputImage('out'))
msks.SetParameterString('exp', '{' + ';'.join(['im1b%d == ' % i + str(ndv) for i in range(1, nd+1)]) + '}')
msks.Execute()
vmsk = otb.Registry.CreateApplication('BandMath')
vmsk.AddImageToParameterInputImageList('il', msks.GetParameterOutputImage('out'))
vmsk.SetParameterString('exp', '+'.join(['im1b%d' % i for i in range(1,nd+1)]) + '!=' + str(nd))
vmsk.Execute()
cmsk = otb.Registry.CreateApplication('BandMathX')
cmsk.AddImageToParameterInputImageList('il', msks.GetParameterOutputImage('out'))
cmsk.AddImageToParameterInputImageList('il', vmsk.GetParameterOutputImage('out'))
cmsk.SetParameterString('exp', '{' + ';'.join(['im1b%d*im2b1' % i for i in range(1,nd+1)]) + '}')
cmsk.Execute()
gf = otb.Registry.CreateApplication('ImageTimeSeriesGapFilling')
gf.SetParameterInputImage('in', sg.GetParameterOutputImage('out'))
gf.SetParameterInputImage('mask', cmsk.GetParameterOutputImage('out'))
gf.SetParameterInt('comp', 1)
gf.SetParameterString('it', 'linear')
gf.SetParameterString('out', output_file)
gf.SetParameterOutputImagePixelType('out', otb.ImagePixelType_int16)
gf.ExecuteAndWriteOutput()
ds = gdal.Open(output_file, gdal.GA_Update)
ds.GetRasterBand(1).SetNoDataValue(ndv)
ds = None
return
def landscape_stratification_metric(fn, output_file, md=[1, 1], duration=365, cbegin=2, cend=5, date_file=None):
if date_file is None:
date_file = os.path.splitext(fn)[0] + '_dates.txt'
ds = gdal.Open(fn)
ndv = ds.GetRasterBand(1).GetNoDataValue()
ds = None
msk = otb.Registry.CreateApplication('ManageNoData')
msk.SetParameterString('in', fn)
msk.SetParameterString('mode', 'buildmask')
msk.Execute()
periods = get_period_intervals(date_file, md, duration)
[print(p) for p in periods]
mts = otb.Registry.CreateApplication('BandMathX')
mts.SetParameterStringList('il', [fn])
mts.SetParameterString('exp', get_mean_expression(periods))
mts.Execute()
pca = otb.Registry.CreateApplication('LandscapeStratificationMetric')
pca.AddImageToParameterInputImageList('ndvits', mts.GetParameterOutputImage('out'))
pca.SetParameterInt('cbegin', cbegin)
pca.SetParameterInt('cend', cend)
pca.SetParameterFloat('bv', ndv)
pca.Execute()
dyn = otb.Registry.CreateApplication('DynamicConvert')
dyn.SetParameterInputImage('in', pca.GetParameterOutputImage('out'))
dyn.SetParameterInt('outmin', 0)
dyn.SetParameterInt('outmax', 2048)
dyn.SetParameterInputImage('mask', msk.GetParameterOutputImage('out'))
dyn.SetParameterOutputImagePixelType('out', otb.ImagePixelType_int16)
dyn.Execute()
fin = otb.Registry.CreateApplication('ManageNoData')
fin.SetParameterInputImage('in', dyn.GetParameterOutputImage('out'))
fin.SetParameterString('mode', 'apply')
fin.SetParameterInputImage('mode.apply.mask', msk.GetParameterOutputImage('out'))
fin.SetParameterString('out', output_file)
fin.SetParameterOutputImagePixelType('out', otb.ImagePixelType_int16)
fin.ExecuteAndWriteOutput()
return
def compute_stratification(metric_fn, th, output_file, cw = 0.5, sw = 0.5):
ds = gdal.Open(metric_fn)
ndv = ds.GetRasterBand(1).GetNoDataValue()
ds = None
seg = otb.Registry.CreateApplication('LSGRM')
seg.SetParameterString('in', metric_fn)
seg.SetParameterString('criterion', 'bs')
seg.SetParameterFloat('criterion.bs.cw', cw)
seg.SetParameterFloat('criterion.bs.sw', sw)
seg.SetParameterFloat('threshold', th)
seg.Execute()
ndt = otb.Registry.CreateApplication('BandMath')
ndt.SetParameterStringList('il', [metric_fn])
ndt.AddImageToParameterInputImageList('il', seg.GetParameterOutputImage('out'))
ndt.SetParameterString('exp', 'im1b1 != ' + str(ndv) + ' ? im2b1 : 0')
ndt.Execute()
fin = otb.Registry.CreateApplication('SimpleVectorization')
fin.SetParameterInputImage('in', ndt.GetParameterOutputImage('out'))
fin.SetParameterString('out', output_file)
fin.ExecuteAndWriteOutput()
return
Markdown is supported
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