Commit c22fd522 authored by Gaetano Raffaele's avatar Gaetano Raffaele
Browse files

WIP: complete refactoring of stratification application.

parent 335fc245
......@@ -22,54 +22,9 @@ from qgis.core import (QgsProcessingAlgorithm,
QgsProcessingParameterRasterDestination)
from qgis import processing
import os
import gdal
import datetime
from shutil import rmtree
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]
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) - delta)
if e is not None:
periods.append((S + s, S + e + 1))
S += s + 1
dates = dates[s + 1:]
return periods
def get_mean_expression(periods):
l = periods[0][1] - periods[0][0]
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) + '}'
from utils import *
def NDVI_subdates_maker(TS_file,datefile,output_folder,begin_date,end_date,context,feedback):
raster_band_list = []
......@@ -93,20 +48,6 @@ def NDVI_subdates_maker(TS_file,datefile,output_folder,begin_date,end_date,conte
f.write(str(date)+'\n')
return out, os.path.join(output_folder,os.path.basename(datefile))
def setNoDataValue(fn, val=None):
ds = gdal.Open(fn, gdal.GA_Update)
if val is None:
ds.GetRasterBand(1).DeleteNoDataValue()
else:
ds.GetRasterBand(1).SetNoDataValue(val)
ds = None
def getNoDataValue(fn):
ds = gdal.Open(fn)
ndv = ds.GetRasterBand(1).GetNoDataValue()
ds = None
return ndv
def moy_NDVI_TS(TS,dates, month, day, duration):
periods = get_period_intervals(dates, [month,day], duration)
return get_mean_expression(periods)
......@@ -337,5 +278,4 @@ class LandscapeStratification(QgsProcessingAlgorithm):
if len(os.listdir(tmp_folder)) == 0:
os.rmdir(tmp_folder)
return {self.OUTPUT: out_layer}
# -*- coding: utf-8 -*-
"""
***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************
"""
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
QgsFeatureSink,
QgsRasterLayer,
QgsProcessingException,
QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSource,
QgsProcessingParameterFolderDestination,
QgsProcessingParameterFileDestination,
QgsProcessingParameterString,
QgsProcessingParameterNumber,
QgsProcessingParameterVectorLayer,
QgsProcessingParameterEnum,
QgsProcessingParameterFeatureSink)
from qgis import processing
import gdal
import glob
import os
import sys
import calendar
import datetime
import shutil
def ModisToYYYYMMDD(raster_list, output):
dates=[]
for path in raster_list :
modisDate = path.split('_')[-2][3:]
year = int(modisDate[:-3])
day = int(modisDate[-3:])
month = 1
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
day = day - calendar.monthrange(year,month)[1]
month = month + 1
month= '0'+str(month) if month < 10 else str(month)
day= '0'+str(day) if day < 10 else str(day)
dates.append(str(year)+str(month)+str(day))
dates.sort()
with open(output,'w') as f:
for date in dates:
f.write(date+'\n')
def setNoDataValue(fn, val=0):
""" setNoDataValue(fn,val = 0)
Sets a value for nodata pixels in a raster GDAL dataset.
INPUT : fn - GDAL dataset file
val - nodata value (def. 0)
OUTPUT : None
"""
ds = gdal.Open(fn, gdal.GA_Update)
for i in range(0, ds.RasterCount):
band = ds.GetRasterBand(i + 1)
if val is None:
band.DeleteNoDataValue()
else:
band.SetNoDataValue(val)
band = None
ds.FlushCache()
ds = None
class TemporalSmoothing(QgsProcessingAlgorithm):
"""
This is an example algorithm that takes a vector layer and
creates a new identical one.
It is meant to be used as an example of how to create your own
algorithms and explain methods and variables used to do it. An
algorithm like this will be available in all elements, and there
is not need for additional work.
All Processing algorithms should extend the QgsProcessingAlgorithm
class.
"""
# Constants used to refer to parameters and outputs. They will be
# used when calling the algorithm from another algorithm, or when
# calling from the QGIS console.
INPUT = 'INPUT'
OUTPUT = 'OUTPUT'
TSINT = 'TSINT'
TSPARAMS = 'TSPARAMS'
def tr(self, string):
"""
Returns a translatable string with the self.tr() function.
"""
return QCoreApplication.translate('Processing', string)
def createInstance(self):
return TemporalSmoothing()
def name(self):
"""
Returns the algorithm name, used for identifying the algorithm. This
string should be fixed for the algorithm, and must not be localised.
The name should be unique within each provider. Names should contain
lowercase alphanumeric characters only and no spaces or other
formatting characters.
"""
return 'TemporalSmoothing'
def displayName(self):
"""
Returns the translated algorithm name, which should be used for any
user-visible display of the algorithm name.
"""
return self.tr('TemporalSmoothing')
def group(self):
"""
Returns the name of the group this algorithm belongs to. This string
should be localised.
"""
return self.tr('TSNDVI')
def groupId(self):
"""
Returns the unique ID of the group this algorithm belongs to. This
string should be fixed for the algorithm, and must not be localised.
The group id should be unique within each provider. Group id should
contain lowercase alphanumeric characters only and no spaces or other
formatting characters.
"""
return 'TSNDVI'
def shortHelpString(self):
"""
Returns a localised short helper string for the algorithm. This string
should provide a basic description about what the algorithm does and the
parameters and outputs associated with it..
"""
return self.tr("Compute the trend of the vegetation index time series")
def initAlgorithm(self, config=None):
"""
Here we define the inputs and output of the algorithm, along
with some other properties.
"""
# We add the input vector features source. It can have any kind of
# geometry.
self.addParameter(
QgsProcessingParameterFolderDestination(
self.INPUT,
self.tr('Input folder')
)
)
self.addParameter(
QgsProcessingParameterEnum(
self.TSINT,
self.tr('Temporal smoothing method : sg/irsg'),
options=['sg','irsg'],
defaultValue='sg'
)
)
self.addParameter(
QgsProcessingParameterString(
self.TSPARAMS,
self.tr('Temporal smoothing parameters : sg.deg sg.rad / irsg.ltdeg irsg.ltrad irsg.stdeg irsg.strad irsg.threshval irsg.threshtime'),
defaultValue='2 2'
))
self.addParameter(
QgsProcessingParameterFileDestination(
self.OUTPUT,
self.tr('Output file')
)
)
def processAlgorithm(self, parameters, context, feedback):
"""
Here is where the processing itself takes place.
"""
TS_folder = self.parameterAsString(
parameters,
self.INPUT,
context
)
output = self.parameterAsString(
parameters,
self.OUTPUT,
context
)
deg=None
rad=None
ltdeg = None
ltrad = None
stdeg = None
strad = None
threshval = None
threshtime = None
options=['sg','irsg']
method = options[parameters['TSINT']]
nb_param = 2 if method == 'sg' else 6
if len(parameters['TSPARAMS'].split(' '))==nb_param :
ts_params = parameters['TSPARAMS'].split(' ')
if len(ts_params) == 2:
deg = ts_params[0]
rad = ts_params[1]
elif len(ts_params) == 6:
ltdeg = ts_params[0]
ltrad = ts_params[1]
stdeg = ts_params[2]
strad = ts_params[3]
threshval = ts_params[4]
threshtime = ts_params[5]
else :
sys.exit('TSPARAMS should have {} parameters with {} smoothing method'.format(nb_param, method))
datefile = os.path.join(os.path.dirname(output), os.path.splitext(parameters['OUTPUT'])[0]+'_dates.txt')
NDVI_TS_files = glob.glob(os.path.join(TS_folder, '*NDVI*.tif'))
NDVI_TS_files.sort()
ModisToYYYYMMDD(NDVI_TS_files,datefile)
VRTOptions = gdal.BuildVRTOptions(separate=True)
if os.path.exists(os.path.join(os.path.dirname(output),'temp')) == False :
os.mkdir(os.path.join(os.path.dirname(output),'temp'))
VRT_file = os.path.join(os.path.dirname(output),'temp','input_ndvi.vrt')
gdal.BuildVRT(VRT_file,NDVI_TS_files,options=VRTOptions)
Smooth_TS = parameters['OUTPUT']
Smooth_app_parameters = {"il": [VRT_file], "dates": datefile, "out": Smooth_TS, "outputpixeltype":2}
#Smooth_app_parameters = {"il": NDVI_TS_files, "dates": datefile, "out": Smooth_TS, "outputpixeltype":2}
if deg != None:
Smooth_app_parameters['interp.sg.deg'] = deg
if rad != None:
Smooth_app_parameters['interp.sg.rad'] = rad
if ltdeg != None:
Smooth_app_parameters['interp.irsg.ltdeg'] = ltdeg
if ltrad != None:
Smooth_app_parameters['interp.irsg.ltrad'] = ltrad
if stdeg != None:
Smooth_app_parameters['interp.irsg.stdeg'] = stdeg
if strad != None:
Smooth_app_parameters['interp.irsg.strad'] = strad
if threshval != None:
Smooth_app_parameters['interp.irsg.threshval'] = threshval
if threshtime != None:
Smooth_app_parameters['interp.irsg.threshtime'] = threshtime
processing.run("otb:TemporalSmoothing", Smooth_app_parameters, context=context, feedback=feedback)
setNoDataValue(Smooth_TS,-3000.0)
shutil.rmtree(os.path.join(os.path.dirname(output),'temp'))
return {'OUTPUT':Smooth_TS,'datefile':datefile}
# -*- coding: utf-8 -*-
"""
***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************
"""
import warnings
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessingAlgorithm,
QgsProcessingParameterFile,
QgsProcessingParameterRasterDestination,
QgsProcessingParameterString,
QgsProcessingParameterEnum)
from qgis import processing
import glob
import os
import sys
import shutil
import datetime
from utils import setNoDataValue, getNoDataValue
class TemporalSmoothing(QgsProcessingAlgorithm):
INPUT = 'INPUT'
OUTPUT = 'OUTPUT'
TSINT = 'TSINT'
TSPARAMS = 'TSPARAMS'
def tr(self, string):
return QCoreApplication.translate('Processing', string)
def createInstance(self):
return TemporalSmoothing()
def name(self):
return 'TemporalSmoothing'
def displayName(self):
return self.tr('TemporalSmoothing')
def group(self):
return self.tr('TSNDVI')
def groupId(self):
return 'TSNDVI'
def shortHelpString(self):
return self.tr("Compute the trend of the vegetation index time series")
def initAlgorithm(self, config=None):
self.addParameter(
QgsProcessingParameterFile(
self.INPUT,
self.tr('Input folder')
)
)
self.addParameter(
QgsProcessingParameterEnum(
self.TSINT,
self.tr('Temporal smoothing method : sg/irsg'),
options=['sg','irsg'],
defaultValue='sg'
)
)
self.addParameter(
QgsProcessingParameterString(
self.TSPARAMS,
self.tr('Temporal smoothing parameters : sg.deg sg.rad / irsg.ltdeg irsg.ltrad irsg.stdeg irsg.strad irsg.threshval irsg.threshtime'),
defaultValue='2 2'
))
self.addParameter(
QgsProcessingParameterRasterDestination(
self.OUTPUT,
self.tr('Output smoothed time series')
)
)
def processAlgorithm(self, parameters, context, feedback):
"""
Here is where the processing itself takes place.
"""
TS_folder = self.parameterAsString(
parameters,
self.INPUT,
context
)
output = self.parameterAsOutputLayer(
parameters,
self.OUTPUT,
context
)
deg=None
rad=None
ltdeg = None
ltrad = None
stdeg = None
strad = None
threshval = None
threshtime = None
options=['sg','irsg']
method = options[parameters['TSINT']]
nb_param = 2 if method == 'sg' else 6
if len(parameters['TSPARAMS'].split(' '))==nb_param :
ts_params = parameters['TSPARAMS'].split(' ')
if len(ts_params) == 2:
deg = ts_params[0]
rad = ts_params[1]
elif len(ts_params) == 6:
ltdeg = ts_params[0]
ltrad = ts_params[1]
stdeg = ts_params[2]
strad = ts_params[3]
threshval = ts_params[4]
threshtime = ts_params[5]
else :
sys.exit('TSPARAMS should have {} parameters with {} smoothing method'.format(nb_param, method))
NDVI_TS_files = sorted(glob.glob(os.path.join(TS_folder, '*NDVI*.tif')))
dates = [x.split('_')[-2][3:] for x in NDVI_TS_files]
sdates = [datetime.datetime.strptime(d, '%Y%j').date().strftime('%Y%m%d') for d in dates]
datefile = os.path.splitext(output)[0] + '_dates.txt'
with open(datefile, 'w') as df:
df.write('\n'.join(sdates) + '\n')
ndv = getNoDataValue(NDVI_TS_files[0])
nd = len(NDVI_TS_files)
'''
VRTOptions = gdal.BuildVRTOptions(separate=True)
VRT_file = os.path.join(os.path.dirname(output),'temp','input_ndvi.vrt')
gdal.BuildVRT(VRT_file,NDVI_TS_files,options=VRTOptions)
'''
tmp_folder = os.path.join(os.path.dirname(output), 'temp')
if not os.path.exists(tmp_folder):
os.mkdir(tmp_folder)
Smooth_TS_pre = os.path.join(tmp_folder, os.path.splitext(os.path.basename(output))[0] + '_tmp.tif')
Smooth_app_parameters = {"il": NDVI_TS_files, "dates": datefile, "out": Smooth_TS_pre, "outputpixeltype":2}
if deg != None:
Smooth_app_parameters['interp.sg.deg'] = deg
if rad != None:
Smooth_app_parameters['interp.sg.rad'] = rad
if ltdeg != None:
Smooth_app_parameters['interp.irsg.ltdeg'] = ltdeg
if ltrad != None:
Smooth_app_parameters['interp.irsg.ltrad'] = ltrad
if stdeg != None:
Smooth_app_parameters['interp.irsg.stdeg'] = stdeg
if strad != None:
Smooth_app_parameters['interp.irsg.strad'] = strad
if threshval != None:
Smooth_app_parameters['interp.irsg.threshval'] = threshval
if threshtime != None:
Smooth_app_parameters['interp.irsg.threshtime'] = threshtime
processing.run("otb:TemporalSmoothing", Smooth_app_parameters, context=context, feedback=feedback)
setNoDataValue(Smooth_TS_pre, ndv)
msks = os.path.join(tmp_folder, os.path.splitext(os.path.basename(output))[0] + '_msks.tif')
msks_exp = '{' + ';'.join(['im1b%d == ' % i + str(ndv) for i in range(1, nd+1)]) + '}'
BMX_params = {'il': Smooth_TS_pre, 'exp': msks_exp, 'out': msks, "outputpixeltype": 0}
processing.run('otb:BandMathX', BMX_params, context=context, feedback=feedback)
vmsk = os.path.join(tmp_folder, os.path.splitext(os.path.basename(output))[0] + '_vmsk.tif')
vmsk_exp = '+'.join(['im1b%d' % i for i in range(1,nd+1)]) + '!=' + str(nd)
BM_params = {'il': [msks], 'exp': vmsk_exp, 'out': vmsk, "outputpixeltype": 0}
processing.run('otb:BandMath', BM_params, context=context, feedback=feedback)
cmsk = os.path.join(tmp_folder, os.path.splitext(os.path.basename(output))[0] + '_cmsk.tif')
cmsk_exp = '{' + ';'.join(['im1b%d*im2b1' % i for i in range(1,nd+1)]) + '}'
BMX_params = {'il': [msks, vmsk], 'exp': cmsk_exp, 'out': cmsk, "outputpixeltype": 0}
processing.run('otb:BandMathX', BMX_params, context=context, feedback=feedback)
gf = output
gf_params = {'in': Smooth_TS_pre, 'out': gf, 'mask': cmsk, 'comp': 1, 'it': 'linear', 'id': datefile, 'od': datefile}
out = processing.run('otb:ImageTimeSeriesGapFilling', gf_params, context=context, feedback=feedback)
setNoDataValue(gf, ndv)
try:
shutil.rmtree(tmp_folder)
except:
warnings.warn('Could not delete temporary folder')
return {self.OUTPUT: out}
import gdal
import datetime
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]
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) - delta)
if e is not None:
periods.append((S + s, S + e + 1))
S += s + 1
dates = dates[s + 1:]
return periods
def get_mean_expression(periods):
l = periods[0][1] - periods[0][0]
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 setNoDataValue(fn, val=None):
ds = gdal.Open(fn, gdal.GA_Update)
if val is None:
ds.GetRasterBand(1).DeleteNoDataValue()
else:
ds.GetRasterBand(1).SetNoDataValue(val)
ds = None
def getNoDataValue(fn):
ds = gdal.Open(fn)
ndv = ds.GetRasterBand(1).GetNoDataValue()
ds = None
return ndv
\ No newline at end of file
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