Commit 3cceec86 authored by Gaetano Raffaele's avatar Gaetano Raffaele
Browse files

ENH: major modifications and corrections.

parent 3aebf146
......@@ -12,22 +12,15 @@
"""
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
QgsFeatureSink,
QgsRasterLayer,
QgsProcessingException,
QgsProcessingAlgorithm,
QgsProcessingParameterFeatureSource,
from qgis.core import (QgsProcessingAlgorithm,
QgsProcessingParameterFolderDestination,
QgsProcessingParameterFile,
QgsProcessingParameterString,
QgsProcessingParameterNumber,
QgsProcessingParameterBoolean,
QgsProcessingParameterVectorLayer,
QgsProcessingParameterEnum,
QgsProcessingParameterFeatureSink)
QgsProcessingContext)
from qgis import processing
import glob
from qgis.utils import iface
import ogr
import os
import sys
......@@ -35,10 +28,6 @@ import gdal
import calendar
import platform
import datetime
import subprocess
import random
import string
import shutil
def find_next_date(datelist, date):
for i,l in enumerate(datelist):
......@@ -153,6 +142,20 @@ def setNoDataValue(fn, val=None):
ds.FlushCache()
ds = None
def getNoDataValue(fn):
""" getNoDataValue(fn,val = 0)
Gets the value for nodata pixels in a raster GDAL dataset.
INPUT : fn - GDAL dataset file
OUTPUT : no-data value
"""
ds = gdal.Open(fn)
band = ds.GetRasterBand(1)
ndv = band.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)
......@@ -260,7 +263,8 @@ class LandscapeStratification(QgsProcessingAlgorithm):
self.addParameter(
QgsProcessingParameterFile(
self.INPUTDATE,
self.tr('Input NDVI dates file')
self.tr('Input NDVI dates file'),
optional=True
)
)
self.addParameter(
......@@ -310,7 +314,9 @@ class LandscapeStratification(QgsProcessingAlgorithm):
self.addParameter(
QgsProcessingParameterVectorLayer(
self.CLIP,
self.tr('Clip vector layer')
self.tr('Clip vector layer'),
optional=True,
defaultValue=None
)
)
......@@ -333,7 +339,7 @@ class LandscapeStratification(QgsProcessingAlgorithm):
self.CEND,
self.tr('Final PCA component (0 for last)'),
type=0,
defaultValue=0))
defaultValue=5))
self.addParameter(
QgsProcessingParameterNumber(
self.SW,
......@@ -413,65 +419,106 @@ class LandscapeStratification(QgsProcessingAlgorithm):
sw=parameters['SW']
cbegin=parameters['CBEGIN']
cend=parameters['CEND']
tmp_name = prefix + 'metrics'
tmp_folder = os.path.join(output_folder,tmp_name)
if not os.path.exists(tmp_folder):
os.mkdir(tmp_folder)
print('Check Datefile')
print(datefile)
if datefile == '':
datefile = os.path.splitext(Smooth_TS)[0] + '_dates.txt'
print(datefile)
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
ndv = getNoDataValue(TS_file)
# Compute mean time series over the period
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"))
maskc = os.path.join(tmp_folder, prefix+"mask_clip.tif")
processing.run("otb:Rasterization", {"in":clip,"im":Smooth_TS,"mode.binary.foreground":1,"out":maskc, "outputpixeltype":0}, context=context, feedback=feedback)
vald = os.path.join(tmp_folder, prefix+"valid.tif")
ND_app_parameters = {"in": TS_file, "out":vald, "mode":"buildmask", 'mode.apply.mask':maskc, "outputpixeltype":0}
processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
mask = os.path.join(tmp_folder, prefix+"mask.tif")
BMX_parameters = {'il': [maskc, vald], 'exp':'im1b1*im2b1', 'out':mask, 'outputpixeltype':0, '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"))
LS_Strat = os.path.join(tmp_folder, prefix+"LandStrat_metric.tif")
LS_app_parameters={"ndvits":Moy_TS, "out": LS_Strat, "outputpixeltype":5}
if ndv is not None:
setNoDataValue(Moy_TS, ndv)
# manage mask (nodata from input U clip mask)
if clip is not None:
maskc = os.path.join(tmp_folder, prefix+"mask_clip.tif")
processing.run("otb:Rasterization", {"in":clip,"im":Smooth_TS,"mode.binary.foreground":1,"out":maskc, "outputpixeltype":0}, context=context, feedback=feedback)
vald = os.path.join(tmp_folder, prefix+"valid.tif")
ND_app_parameters = {"in": TS_file, "out":vald, "mode":"buildmask", 'mode.apply.mask':maskc, "outputpixeltype":0}
processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
mask = os.path.join(tmp_folder, prefix+"mask.tif")
BMX_parameters = {'il': [maskc, vald], 'exp':'im1b1*im2b1', 'out':mask, 'outputpixeltype':0, '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"))
os.remove(maskc)
os.remove(vald)
else:
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)
# 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
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)
LS_Strat = os.path.join(tmp_folder, prefix + "LandStrat_metric.tif")
ND_app_parameters = {'in': LS_Strat_pre, 'out': LS_Strat, 'mode': 'apply', 'mode.apply.mask': mask,
'outputpixeltype': 5}
processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
setNoDataValue(LS_Strat, ndv)
LS_Strat_norm = os.path.join(tmp_folder, prefix+"LandStrat_metric_norm.tif")
DC_app_parameters={"in":LS_Strat, "out":LS_Strat_norm, "outmin":0, "outmax":2048, "outputpixeltype":2}
DC_app_parameters={"in":LS_Strat, "out":LS_Strat_norm, "mask":mask, "outmin":0, "outmax":2048, "outputpixeltype":2}
processing.run("otb:DynamicConvert", DC_app_parameters, context=context, feedback=feedback)
#setNoDataValue(LS_Strat_norm,-9999)
NDoutput = os.path.join(tmp_folder,prefix+"LandStrat_metric_norm_masked.tif")
ND_app_parameters = {'in': LS_Strat_norm, 'out':NDoutput, 'mode':'apply', 'mode.apply.mask':mask, 'outputpixeltype':2}
LS_Strat = os.path.join(tmp_folder,prefix+"LandStrat_metric_norm_masked.tif")
ND_app_parameters = {'in': LS_Strat_norm, 'out':LS_Strat, 'mode':'apply', 'mode.apply.mask':mask, 'outputpixeltype':2}
processing.run('otb:ManageNoData', ND_app_parameters, context=context, feedback=feedback)
LS_Strat = NDoutput
setNoDataValue(LS_Strat, ndv)
os.remove(LS_Strat_pre)
os.remove(LS_Strat_norm)
GRM_output = os.path.join(tmp_folder,prefix+"LandStrat_map.tif")
GRM_app_parameters = {'in':LS_Strat, 'threshold':float(threshold), 'criterion':'bs', 'out':GRM_output}
GRM_output_pre = os.path.join(tmp_folder,prefix+"LandStrat_map_pre.tif")
GRM_app_parameters = {'in':LS_Strat, 'threshold':float(threshold), 'criterion':'bs', 'out':GRM_output_pre}
if cw != 0.5:
GRM_app_parameters['criterion.bs.cw']= cw
if sw != 0.5:
GRM_app_parameters['criterion.bs.sw']= sw
processing.run("otb:LSGRM", GRM_app_parameters, context=context, feedback=feedback)
GRM_output = os.path.join(tmp_folder, prefix + "LandStrat_map.tif")
BMX_parameters = {'il': [GRM_output_pre, mask], 'exp': 'im1b1*im2b1', 'out': GRM_output, 'outputpixeltype': 3,
'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"))
os.remove(GRM_output_pre)
setNoDataValue(GRM_output, 0)
GRM_vec_output = os.path.join(output_folder,prefix+"LandscapeStratification.shp")
if os.path.exists(GRM_vec_output):
drv = ogr.GetDriverByName('ESRI Shapefile')
drv.DeleteDataSource(GRM_vec_output)
processing.run("gdal:polygonize", {'INPUT':GRM_output, 'BAND':1, 'OUTPUT':GRM_vec_output}, context=context, feedback=feedback)
nm = os.path.basename(os.path.splitext(GRM_vec_output)[0])
context.addLayerToLoadOnCompletion(GRM_vec_output,
QgsProcessingContext.LayerDetails(name=nm, project=context.project()))
return {'OUT':GRM_vec_output}
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