LandscapeStratificationMetric.py 11 KB
Newer Older
Gaetano Raffaele's avatar
Gaetano Raffaele committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# -*- 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
15
from qgis.core import (QgsProcessingAlgorithm,
Gaetano Raffaele's avatar
Gaetano Raffaele committed
16
17
18
                       QgsProcessingParameterFile,
                       QgsProcessingParameterString,
                       QgsProcessingParameterNumber,
19
20
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterRasterDestination)
Gaetano Raffaele's avatar
Gaetano Raffaele committed
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
from qgis import processing
import os
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)

40
def get_period_intervals(date_file, md=[1, 1], duration=365):
Gaetano Raffaele's avatar
Gaetano Raffaele committed
41
42
43
44
45
46
47
    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
48
    delta = dates[-1] - dates[-2]
Gaetano Raffaele's avatar
Gaetano Raffaele committed
49
    periods = []
50
    s, e = 0, 0
Gaetano Raffaele's avatar
Gaetano Raffaele committed
51
52
53
54
    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:
55
            e = find_next_date(dates, dates[s] + datetime.timedelta(days=duration) - delta)
Gaetano Raffaele's avatar
Gaetano Raffaele committed
56
            if e is not None:
57
58
59
                periods.append((S + s, S + e + 1))
                S += s + 1
                dates = dates[s + 1:]
Gaetano Raffaele's avatar
Gaetano Raffaele committed
60
    return periods
61

Gaetano Raffaele's avatar
Gaetano Raffaele committed
62
def get_mean_expression(periods):
63
    l = periods[0][1] - periods[0][0]
Gaetano Raffaele's avatar
Gaetano Raffaele committed
64
    expr = []
65
66
    for i in range(1, l + 1):
        s = [p[0] + i for p in periods]
Gaetano Raffaele's avatar
Gaetano Raffaele committed
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
        expr.append('(' + '+'.join(['im1b%d' % x for x in s]) + ')/' + str(len(s)))

    return '{' + ';'.join(expr) + '}'

def NDVI_subdates_maker(TS_file,datefile,output_folder,begin_date,end_date,context,feedback):
    raster_band_list = []
    dates_list= []
    with open(datefile,'r') as f:
        for i,date in enumerate(f):
            if int(date)>=begin_date and int(date) <=end_date:
                raster_band_list.append(i+1)
                dates_list.append(int(date))
    exp=''
    for i in raster_band_list :
        exp+='im1b'+str(i)+';'
    exp=exp[:-1]
    out = os.path.join(output_folder,"sub_ndvi_ts.tif")
    BMX_parameters = {'il':[TS_file],'exp':exp, 'out':out, 'outcontext':os.path.join(output_folder,"outcontext.txt")}
    processing.run('otb:BandMathX',BMX_parameters, context=context, feedback=feedback)
    os.remove(os.path.join(output_folder,"outcontext.txt"))
    
    with open(os.path.join(output_folder,os.path.basename(datefile)),'w') as f:
        for date in dates_list:
            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)
95
96
97
98
    if val is None:
        ds.GetRasterBand(1).DeleteNoDataValue()
    else:
        ds.GetRasterBand(1).SetNoDataValue(val)
Gaetano Raffaele's avatar
Gaetano Raffaele committed
99
100
    ds = None

101
102
def getNoDataValue(fn):
    ds = gdal.Open(fn)
103
    ndv = ds.GetRasterBand(1).GetNoDataValue()
104
105
106
    ds = None
    return ndv

Gaetano Raffaele's avatar
Gaetano Raffaele committed
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 moy_NDVI_TS(TS,dates, month, day, duration):
    periods = get_period_intervals(dates, [month,day], duration)
    return get_mean_expression(periods)

class LandscapeStratification(QgsProcessingAlgorithm):

    INPUT = 'INPUT'
    INPUTDATE = 'INPUTDATE'
    BEGDATE = 'BEGDATE'
    ENDDATE = 'ENDDATE'
    STARTMONTH = 'STARTMONTH'
    STARTDAY = 'STARTDAY'
    DURATION = 'DURATION'
    OUTPUT = 'OUTPUT'
    PREFIX = 'PREFIX'
    CBEGIN = 'CBEGIN'
    CEND = 'CEND'

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return LandscapeStratification()

    def name(self):
        return 'LandscapeStratification'

    def displayName(self):
        return self.tr('LandscapeStratification')

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

    def groupId(self):
        return 'TSNDVI'

    def shortHelpString(self):
        return self.tr("Compute a landscape stratification metric")

    def initAlgorithm(self, config=None):
        self.addParameter(
148
            QgsProcessingParameterRasterLayer(
Gaetano Raffaele's avatar
Gaetano Raffaele committed
149
150
151
152
153
154
155
156
                self.INPUT,
                self.tr('Input NDVI stack')
            )
        )

        self.addParameter(
            QgsProcessingParameterFile(
                self.INPUTDATE,
157
158
                self.tr('Input NDVI dates file'),
                optional=True
Gaetano Raffaele's avatar
Gaetano Raffaele committed
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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
            )
        )
        self.addParameter(
            QgsProcessingParameterString(
                self.BEGDATE,
                self.tr('Beginning date for selection (YYYYMMDD)'),
                defaultValue='20000101'
            )
        )
        self.addParameter(
            QgsProcessingParameterString(
                self.ENDDATE,
                self.tr('End date for selection (YYYYMMDD)'),
                defaultValue=datetime.date.today().strftime("%Y%m%d")
            )
        )
        self.addParameter(
            QgsProcessingParameterNumber(
                self.STARTMONTH,
                self.tr('Starting month for multi-annual mean'),
                type=0,
                defaultValue=1,
                minValue=1,
                maxValue=12
            )
        )
        self.addParameter(
            QgsProcessingParameterNumber(
                self.STARTDAY,
                self.tr('Starting day for multi-annual mean'),
                type=0,
                defaultValue=1,
                minValue=1,
                maxValue=31
            )
        )
        self.addParameter(
            QgsProcessingParameterNumber(
                self.DURATION,
                self.tr('Period length for multi-annual mean'),
                type=0,
                defaultValue=365,
                minValue=1
            )
        )

        self.addParameter(
            QgsProcessingParameterNumber(
                self.CBEGIN,
                self.tr('Starting PCA component'),
                type=0,
                defaultValue=2))
211

Gaetano Raffaele's avatar
Gaetano Raffaele committed
212
213
214
215
216
        self.addParameter(
            QgsProcessingParameterNumber(
                self.CEND,
                self.tr('Final PCA component (0 for last)'),
                type=0,
217
                defaultValue=5))
Gaetano Raffaele's avatar
Gaetano Raffaele committed
218
219

        self.addParameter(
220
            QgsProcessingParameterRasterDestination(
Gaetano Raffaele's avatar
Gaetano Raffaele committed
221
                self.OUTPUT,
222
                self.tr('Output landscape stratification metric')
Gaetano Raffaele's avatar
Gaetano Raffaele committed
223
224
225
226
            )
        )

    def processAlgorithm(self, parameters, context, feedback):
227

228
229
230
231
232
        Smooth_TS = self.parameterAsString(
            parameters,
            self.INPUT,
            context
        )
Gaetano Raffaele's avatar
Gaetano Raffaele committed
233
234
235
236
237
238
239
        
        datefile = self.parameterAsString(
            parameters,
            self.INPUTDATE,
            context
        )

240
241
242
243
244
245
        output_raster = self.parameterAsOutputLayer(
            parameters,
            self.OUTPUT,
            context
        )

246
247
        output_folder = os.path.dirname(output_raster)

248
        prefix = os.path.splitext(os.path.basename(output_raster))[0] + '_'
Gaetano Raffaele's avatar
Gaetano Raffaele committed
249
250
        begin_date = parameters['BEGDATE'] if parameters['BEGDATE'] != None else 20000101
        end_date = parameters['ENDDATE'] if parameters['ENDDATE'] != None else int(datetime.date.today().strftime("%Y%m%d"))
251

252
253
254
255
        cbegin = parameters['CBEGIN']
        cend = parameters['CEND']

        tmp_name = prefix + 'temp'
Gaetano Raffaele's avatar
Gaetano Raffaele committed
256
257
258
259
        tmp_folder = os.path.join(output_folder,tmp_name)
        if not os.path.exists(tmp_folder):
            os.mkdir(tmp_folder)

260
261
262
        if datefile == '':
            datefile = os.path.splitext(Smooth_TS)[0] + '_dates.txt'

Gaetano Raffaele's avatar
Gaetano Raffaele committed
263
264
265
266
267
        if int(begin_date)!=20000101 or int(end_date)!=int(datetime.date.today().strftime("%Y%m%d")):
            TS_file, datefile = NDVI_subdates_maker(Smooth_TS,datefile,tmp_folder,int(begin_date),int(end_date),context,feedback)
        else :
            TS_file = Smooth_TS

268
269
270
        ndv = getNoDataValue(TS_file)

        # Compute mean time series over the period
Gaetano Raffaele's avatar
Gaetano Raffaele committed
271
272
273
274
275
        exp = moy_NDVI_TS(TS_file,datefile,parameters['STARTMONTH'], parameters['STARTDAY'], parameters['DURATION'])
        Moy_TS = os.path.join(tmp_folder, prefix+"LandStrat_mean.tif")
        BMX_parameters = {'il':[TS_file],'exp':exp, 'out':Moy_TS, "outputpixeltype":2, 'outcontext':os.path.join(output_folder,"outcontext.txt")}
        processing.run("otb:BandMathX", BMX_parameters, context=context, feedback=feedback)
        os.remove(os.path.join(output_folder,"outcontext.txt"))
276
277
278
        if ndv is not None:
            setNoDataValue(Moy_TS, ndv)

279
280
281
282
        mask = os.path.join(tmp_folder, prefix + "mask.tif")
        ND_app_parameters = {"in": TS_file, "out": mask, "mode": "buildmask", 'mode.apply.mask': TS_file,
                             "outputpixeltype": 0}
        processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
283
284
285
286
287
288

        # Compute metric using no-data value
        LS_Strat_pre = os.path.join(tmp_folder, prefix+"LandStrat_metric_pre.tif")
        LS_app_parameters = {"ndvits": Moy_TS, "out": LS_Strat_pre, "outputpixeltype": 5}
        if ndv is not None:
            LS_app_parameters["bv"] = ndv
Gaetano Raffaele's avatar
Gaetano Raffaele committed
289
290
291
292
293
        if cbegin != None :
            LS_app_parameters["cbegin"]=cbegin
        if cend != None :
            LS_app_parameters["cend"]=cend
        processing.run("otb:LandscapeStratificationMetric", LS_app_parameters, context=context, feedback=feedback)
294

295
296
        LS_Strat_raw = os.path.join(tmp_folder, prefix + "LandStrat_metric_raw.tif")
        ND_app_parameters = {'in': LS_Strat_pre, 'out': LS_Strat_raw, 'mode': 'apply', 'mode.apply.mask': mask,
297
298
                             'outputpixeltype': 5}
        processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
299
        setNoDataValue(LS_Strat_raw, ndv)
Gaetano Raffaele's avatar
Gaetano Raffaele committed
300
301
        
        LS_Strat_norm = os.path.join(tmp_folder, prefix+"LandStrat_metric_norm.tif")
302
        DC_app_parameters={"in":LS_Strat_raw, "out":LS_Strat_norm, "mask":mask, "outmin":0, "outmax":2048, "outputpixeltype":2}
Gaetano Raffaele's avatar
Gaetano Raffaele committed
303
304
        processing.run("otb:DynamicConvert", DC_app_parameters, context=context, feedback=feedback)

305
        LS_Strat = output_raster
306
        ND_app_parameters = {'in': LS_Strat_norm, 'out':LS_Strat, 'mode':'apply', 'mode.apply.mask':mask, 'outputpixeltype':2}
307
        out_layer = processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
308
309
        setNoDataValue(LS_Strat, ndv)

310
        return {self.OUTPUT: out_layer}