TemporalSmoothing.py 7.45 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# -*- 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
Gaetano Raffaele's avatar
Gaetano Raffaele committed
23
import gdal
24
25
26
27
28
import os
import sys
import shutil
import datetime

29
sys.path.append(os.path.dirname(__file__))
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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):
49
        return self.tr('Temporal Smoothing')
50
51
52
53
54
55
56
57

    def group(self):
        return self.tr('TSNDVI')

    def groupId(self):
        return 'TSNDVI'

    def shortHelpString(self):
58
        return self.tr("Apply Savitzky-Golay filtering to a single channel time series (e.g. NDVI).")
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

    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)

        tmp_folder = os.path.join(os.path.dirname(output), 'temp')
        if not os.path.exists(tmp_folder):
            os.mkdir(tmp_folder)
Gaetano Raffaele's avatar
Gaetano Raffaele committed
148
149
150
151
152
            
        ts_vrt_opt = gdal.BuildVRTOptions(separate=True)
        ts_vrt = os.path.join(tmp_folder, os.path.splitext(os.path.basename(TS_folder))[0] + '_list.vrt')
        gdal.BuildVRT(ts_vrt, NDVI_TS_files, options=ts_vrt_opt)
        
153
        Smooth_TS_pre = os.path.join(tmp_folder, os.path.splitext(os.path.basename(output))[0] + '_tmp.tif')
Gaetano Raffaele's avatar
Gaetano Raffaele committed
154
        Smooth_app_parameters = {"il": [ts_vrt], "dates": datefile, "out": Smooth_TS_pre, "outputpixeltype":2}
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
        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
190
        gf_params = {'in': Smooth_TS_pre, 'out': gf, 'mask': cmsk, 'comp': 1, 'it': 'linear', 'id': datefile, 'od': datefile, "outputpixeltype": 2}
191
192
193
194
195
196
197
198
199
        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}