Skip to content

Commit

Permalink
Update agg_average.py
Browse files Browse the repository at this point in the history 
agg_filter finish
  • Loading branch information
@cefect
cefect committed Oct 24, 2022
1 parent cbd2000 commit efd613c
Showing 1 changed file with 80 additions and 27 deletions.
107 changes: 80 additions & 27 deletions floodrescaler/processing/agg_average.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
QgsProcessingParameterString,
QgsRasterLayer ,
QgsCoordinateTransformContext,
QgsMapLayerStore

)

Expand Down Expand Up @@ -74,7 +75,7 @@ def displayName(self):
Returns the translated algorithm name, which should be used for any
user-visible display of the algorithm name.
"""
return self.tr('Aggregation with Averaging')
return self.tr('Aggregation via Averaging')

def group(self):
"""
Expand Down Expand Up @@ -237,9 +238,13 @@ def _agg_simple(self, inp, scale,
'TARGET_RESOLUTION' : scale}

#gives a filepath regardless
return processing.run('gdal:warpreproject', pars_d,
ofp = processing.run('gdal:warpreproject', pars_d,
context=context, feedback=feedback, is_child_algorithm=True)['OUTPUT']

if not os.path.exists(ofp):
raise QgsProcessingException('gdal:warpreproject failed to get a result for \n%s'%pars_d['INTPUT'])

return ofp


def _gdal_calc_1(self, inp, formula, output='TEMPORARY_OUTPUT', **kwargs):
Expand All @@ -253,18 +258,7 @@ def _gdal_calc_1(self, inp, formula, output='TEMPORARY_OUTPUT', **kwargs):
#'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'EXTRA' : '', 'OPTIONS' : '',
}, **kwargs)

def _gdal_calc_add(self, inpA, inpB, output='TEMPORARY_OUTPUT', **kwargs):
"""add two rasters together"""

pars_d = {
'INPUT_A' : inpA, 'BAND_A' : 1, 'INPUT_B' : inpB, 'BAND_B' : 1,
'FORMULA':'A+B',
'NO_DATA' : -9999, 'OUTPUT' : output, 'RTYPE' : 5,
#'BAND_B' : None, 'BAND_C' : None, 'BAND_D' : None, 'BAND_E' : None, 'BAND_F' : None,
#'INPUT_B' : None, 'INPUT_C' : None, 'INPUT_D' : None, 'INPUT_E' : None, 'INPUT_F' : None, 'EXTRA' : '', 'OPTIONS' : '',
}

return self._gdal_calc(pars_d, **kwargs)


def _gdal_calc(self, pars_d, context=None, feedback=None):

Expand All @@ -290,6 +284,7 @@ def get_out(attn):
#=======================================================================
# compute DEM and WSh
#=======================================================================
feedback.pushInfo('computing DEM')
# simple DEM aggregate
dem2_fp = self._agg_simple(dem1, scale,
# output=params['OUTPUT_DEM'],
Expand All @@ -313,7 +308,12 @@ def get_out(attn):
wsh2_maskd_fp = self._gdal_calc_1(wsh2_fp,'A*(A!=0)/(A!=0)', **cf_kwargs)

#DEm + WSH
wse2_fp = self._gdal_calc_add(wsh2_maskd_fp, dem2_fp, output=get_out(self.OUTPUT_WSE), **cf_kwargs)
wse2_fp = self._gdal_calc({
'INPUT_A' : wsh2_maskd_fp, 'BAND_A' : 1, 'INPUT_B' : dem2_fp, 'BAND_B' : 1,
'FORMULA':'A+B',
'NO_DATA' : -9999, 'OUTPUT' : get_out(self.OUTPUT_WSE), 'RTYPE' : 5,
}, **cf_kwargs)




Expand All @@ -330,42 +330,86 @@ def agg_filter(self, params, dem1, wsh1, wse1, scale, context=None, feedback=Non
def get_out(attn):
return self.parameterAsOutputLayer(params, getattr(self, attn), context)

mstore = QgsMapLayerStore()

#=======================================================================
# compute DEM and WSE
#=======================================================================

# simple DEM aggregate
feedback.pushInfo('computing DEM')
dem2_fp = self._agg_simple(dem1, scale, output=get_out(self.OUTPUT_DEM), **cf_kwargs)

#simple WSH aggregate
wse2A_fp = self._agg_simple(wsh1, scale,**cf_kwargs)
feedback.pushInfo('computing raw WSE')
wse2A_fp = self._agg_simple(wse1, scale,**cf_kwargs)

if feedback.isCanceled():
return {}


#=======================================================================
# filter WSE
#=======================================================================
#=======================================================================

feedback.pushInfo('filtering WSE')

"""throwing a permissions error
wse2_mask_fp = self._gdal_calc({
'INPUT_A' : dem2_fp, 'BAND_A' : 1, 'INPUT_B' : wse2A_fp, 'BAND_B' : 1,
#'FORMULA':'B/(A<B)', #WSE / (wet=1, dry=0)
'FORMULA':'B<A', #wet=1, dry=0
'FORMULA':'B>A', #wet=1, dry=0
'NO_DATA' : -9999, 'OUTPUT' :'TEMPORARY_OUTPUT', 'RTYPE' : 5
}, **cf_kwargs)
}, **cf_kwargs)"""

#=======================================================================
# build mask
#=======================================================================
dem2 = QgsRasterLayer(dem2_fp, 'DEM2')
mstore.addMapLayer(dem2)

wse2A = QgsRasterLayer(wse2A_fp, 'WSE2a')
mstore.addMapLayer(wse2A)

with RasterCalc(dem2, feedback=feedback) as rcalc:
dem_rc = rcalc.ref_rcentry
wse_rc = rcalc.add_rcentry(wse2A)

wse2_mask_fp = rcalc.processCalculation('{wse}>{dem}'.format(wse=wse_rc.ref, dem=dem_rc.ref))


#=======================================================================
# apply mask
#=======================================================================
wse2_fp = self._gdal_calc({
'INPUT_A' : wse2_mask_fp, 'BAND_A' : 1, 'INPUT_B' : wse2A_fp, 'BAND_B' : 1,
'FORMULA':'B/A', #WSE / (wet=1, dry=0)
'INPUT_A' : wse2A_fp, 'BAND_A' : 1, 'INPUT_B' : wse2_mask_fp, 'BAND_B' : 1,
'FORMULA':'A/B', #WSE / (wet=1, dry=0)
'NO_DATA' : -9999, 'OUTPUT' : get_out(self.OUTPUT_WSE), 'RTYPE' : 5
}, **cf_kwargs)



if feedback.isCanceled():
return {}
#=======================================================================
# add WSH
# WSE-DEM = WSH
#=======================================================================
wsh2_fp = self._gdal_calc_add(wse2_fp, dem2_fp, output=get_out(self.OUTPUT_WSH), **cf_kwargs)
feedback.pushInfo('computing WSH')
wsh2A_fp = self._gdal_calc({
'INPUT_A' : wse2_fp, 'BAND_A' : 1, 'INPUT_B' : dem2_fp, 'BAND_B' : 1,
'FORMULA':'A-B',
'NO_DATA' : -9999, 'OUTPUT' :'TEMPORARY_OUTPUT', 'RTYPE' : 5,
}, **cf_kwargs)

#fill zeros
wsh2_fp = processing.run('native:fillnodata',
{ 'BAND' : 1,'FILL_VALUE' : 0.0,'INPUT' : wsh2A_fp,'OUTPUT' : get_out(self.OUTPUT_WSH)},
is_child_algorithm=True, **cf_kwargs)['OUTPUT']

#=======================================================================
# wrap
#=======================================================================

mstore.removeAllMapLayers()
return {self.OUTPUT_DEM:dem2_fp, self.OUTPUT_WSH:wsh2_fp, self.OUTPUT_WSE:wse2_fp}


Expand Down Expand Up @@ -415,13 +459,17 @@ def add_rcentry(self, rlay, bandNumber=1):
self.rasterEntries.append(rcentry)
return rcentry

def rcalc(self, formula, output,
def processCalculation(self, formula,
output=None,
ref_lay=None, rasterEntries=None):
#=======================================================================
# defaults
#=======================================================================
if ref_lay is None: ref_lay = self.ref_lay
if rasterEntries is None: rasterEntries = self.rasterEntries
if output is None:
#couldnt figur eout how to access the processing temp dir
output=os.path.join(os.environ['TEMP'], 'processCalculationResult.tif')

assert isinstance(rasterEntries, list)
#=======================================================================
Expand All @@ -433,7 +481,7 @@ def rcalc(self, formula, output,
ofp=output,
outputExtent=ref_lay.extent(),
outputFormat='GTiff',
crs=self.qproj.crs(),
crs=ref_lay.crs(),
nOutputColumns=ref_lay.width(),
nOutputRows=ref_lay.height(),
crsTrnsf=QgsCoordinateTransformContext(),
Expand All @@ -458,13 +506,18 @@ def rcalc(self, formula, output,
if not result == 0:
raise QgsProcessingException('formula=%s failed w/ \n %s' % (formula, rcalc.lastError()))

if not os.path.exists(output):
raise QgsProcessingException(f'failed to get a result w/ \n {formula}')

self.feedback.pushInfo(f'finished on {output}')

return output

def __enter__(self, *args, **kwargs):
return self


def __exit__(self, *args, **kwargs):
return

def assert_extent_equal(left, right, msg='',):
""" extents check"""
Expand Down

0 comments on commit efd613c

Please sign in to comment.