#!/usr/bin/env python # -*- coding: utf-8 -*- ############################################################################ # # MODULE: r.mblend # # AUTHOR(S): Luís Moreira de Sousa # # PURPOSE: Blends two rasters of different spatial resolution. # # COPYRIGHT: (c) 2017 Luís Moreira de Sousa # # This programme is released under the European Union Public # Licence v 1.1. Please consult the LICENCE file for details. # ############################################################################# #%module #% description: Blends two rasters of different spatial resolution. #% keyword: raster #% keyword: resolution #%end #%option HIGH #% key: high #% description: High resolution input raster. #%end #%option LOW #% key: low #% description: Low resolution input raster. #%end #%option OUTPUT #% key: output #% description: Name of output raster. #%end #%option FAR_EDGE #% key: far_edge #% key_desc: value #% type: double #% description: Percentage of distance to high resolution raster used to determine far edge. Number between 0 and 100. #% answer: 95 #% multiple: no #% required: no #%end #%option INTER_POINTS #% key: inter_points #% key_desc: value #% type: integer #% description: Number of points to use in interpolation. A higher number produces a smoother result but requires a lengthier computation. #% answer: 50 #% multiple: no #% required: no #%end #%flag #% key: a #% description: Assign the average difference between the two rasters to the far edge (instead of zero). #%end import os import atexit import math from time import gmtime, strftime import grass.script as gscript index = 0 far_edge_value = '0' d_max = None TMP_MAPS = [] WEIGHT_MAX = 10000 COL_VALUE = 'value' COL_FLAG = 'flag' def getTemporaryIdentifier(): global index global TMP_MAPS id = 'tmp_' + str(os.getpid()) + str(index) index = index + 1 TMP_MAPS.append(id) return id def cleanup(): while len(TMP_MAPS) > 0: gscript.run_command('g.remove', type='all', name=TMP_MAPS.pop(), flags='f', quiet=True) def compute_d_max(region): global d_max print("Region:\n" + str(region)) d_max = math.sqrt(math.pow(region['w'] - region['e'], 2) + math.pow(region['n'] - region['s'], 2)) / 100 def main(): global far_edge_value global d_max options, flags = gscript.parser() high = options['high'] low = options['low'] output = options['output'] far_edge = float(options['far_edge']) inter_points = int(options['inter_points']) use_average_differences = flags['a'] if(high is None or high == ""): gscript.error(_("[r.mblend] ERROR: high is a mandatory parameter.")) exit() if(low is None or low == ""): gscript.error(_("[r.mblend] ERROR: low is a mandatory parameter.")) exit() if(output is None or output == ""): gscript.error(_("[r.mblend] ERROR: output is a mandatory parameter.")) exit() if(far_edge < 0 or far_edge > 100): gscript.error(_("[r.mblend] ERROR: far_edge must be a percentage", " between 0 and 100.")) exit() if(inter_points < 0): gscript.error(_("[r.mblend] ERROR: inter_points must be a positive", " integer.")) exit() # Set the region to the two input rasters gscript.run_command('g.region', raster=high + "," + low) # Determine cell side region = gscript.region() if region['nsres'] > region['ewres']: cell_side = region['nsres'] else: cell_side = region['ewres'] compute_d_max(region) # Make cell size compatible low_res_inter = getTemporaryIdentifier() gscript.message(_("[r.mblend] Resampling low resolution raster to higher" + " resolution")) gscript.run_command('r.resamp.interp', input=low, output=low_res_inter, method='nearest') # Obtain extent to interpolate low_extent_rast = getTemporaryIdentifier() high_extent_rast = getTemporaryIdentifier() low_extent = getTemporaryIdentifier() high_extent = getTemporaryIdentifier() interpol_area = getTemporaryIdentifier() gscript.message(_("[r.mblend] Multiplying low resolution by zero")) gscript.mapcalc(low_extent_rast + ' = ' + low + ' * 0') gscript.message(_("[r.mblend] Multiplying high resolution by zero")) gscript.mapcalc(high_extent_rast + ' = ' + high + ' * 0') gscript.message(_("[r.mblend] Computing extent of low resolution")) gscript.run_command('r.to.vect', input=low_extent_rast, output=low_extent, type='area') gscript.message(_("[r.mblend] Computing extent of high resolution")) gscript.run_command('r.to.vect', input=high_extent_rast, output=high_extent, type='area') gscript.message(_("[r.mblend] Computing area to interpolate")) gscript.run_command('v.overlay', ainput=low_extent, binput=high_extent, output=interpol_area, operator='not') # Compute difference between the two rasters and vectorise to points interpol_area_buff = getTemporaryIdentifier() diff = getTemporaryIdentifier() diff_points_edge = getTemporaryIdentifier() gscript.mapcalc(diff + ' = ' + high + ' - ' + low_res_inter) gscript.message(_("[r.mblend] Computing buffer around interpolation area")) gscript.run_command('v.buffer', input=interpol_area, output=interpol_area_buff, type='area', distance=cell_side) gscript.message(_("[r.mblend] Vectorising differences between input" + " rasters")) gscript.run_command('r.mask', vector=interpol_area_buff) gscript.run_command('r.to.vect', input=diff, output=diff_points_edge, type='point') gscript.run_command('r.mask', flags='r') # Compute average of the differences if flag -a was passed if use_average_differences: p = gscript.pipe_command('r.univar', map=diff) result = {} for line in p.stdout: vector = line.split(": ") if vector[0] == "mean": print("Found it: " + vector[1]) far_edge_value = vector[1] p.wait() # Get points in low resolution farther away from high resolution raster dist_high = getTemporaryIdentifier() weights = getTemporaryIdentifier() interpol_area_points = getTemporaryIdentifier() pre_interpol_area_points = getTemporaryIdentifier() weight_points = getTemporaryIdentifier() interpol_area_in_buff = getTemporaryIdentifier() weight_points_all_edges = getTemporaryIdentifier() weight_points_edge = getTemporaryIdentifier() # 1. Distance to High resolution raster gscript.message(_("[r.mblend] Computing distance to high resolution" + " raster")) gscript.run_command('r.grow.distance', input=high, distance=dist_high) # 2. Rescale to the interval [0,10000]: these are the weights gscript.message(_("[r.mblend] Rescaling distance to [0,10000] interval")) gscript.run_command('r.rescale', input=dist_high, output=weights, to='0,' + str(WEIGHT_MAX)) # 3. Extract points from interpolation area border gscript.message(_("[r.mblend] Extract points from interpolation area " + "boundary")) gscript.run_command('v.to.points', input=interpol_area, output=pre_interpol_area_points, type='boundary', dmax=d_max, layer='-1') gscript.message(_("[r.mblend] Copying features to layer 1")) gscript.run_command('v.category', input=pre_interpol_area_points, output=interpol_area_points, option='chlayer', layer='2,1') gscript.message(_("[r.mblend] Linking attribute table to layer 1")) gscript.run_command('v.db.connect', map=interpol_area_points, table=interpol_area_points, layer='1', flags='o') # 4. Query distances to interpolation area points gscript.message(_("[r.mblend] Querying distances raster")) gscript.run_command('v.what.rast', map=interpol_area_points, raster=weights, column=COL_VALUE) # 5. Select those with higher weights cut_off = str(far_edge / 100 * WEIGHT_MAX) gscript.message(_("[r.mblend] Selecting far edge points (using cut-off" + " percentage)")) gscript.run_command('v.extract', input=interpol_area_points, output=weight_points_edge, where=COL_VALUE + '>' + cut_off) # Merge the two point edges and set low res edge to zero points_edges = getTemporaryIdentifier() gscript.message(_("[r.mblend] Dropping extra column from far edge")) gscript.run_command('v.db.dropcolumn', map=weight_points_edge, layer='1', columns='along') gscript.message(_("[r.mblend] Setting far edge weights to zero")) gscript.run_command('v.db.update', map=weight_points_edge, column=COL_VALUE, value=far_edge_value) gscript.message(_("[r.mblend] Patching the two edges")) gscript.run_command('v.patch', input=diff_points_edge + ',' + weight_points_edge, output=points_edges, flags='e') # Interpolate smoothing raster smoothing = getTemporaryIdentifier() interpol_area_rst = getTemporaryIdentifier() # Consign region to interpolation area gscript.run_command('g.region', vector=interpol_area_buff) gscript.message(_("[r.mblend] Interpolating smoothing surface. This" + " might take a while...")) gscript.run_command('v.surf.idw', input=points_edges, column=COL_VALUE, output=smoothing, power=2, npoints=inter_points) # Reset region to full extent gscript.run_command('g.region', raster=high + "," + low) # Apply stitching smooth_low_res = getTemporaryIdentifier() # Sum to low res gscript.message(_("[r.mblend] Applying smoothing surface")) gscript.mapcalc(smooth_low_res + ' = ' + low_res_inter + ' + ' + smoothing) # Add both rasters try: gscript.message(_("[r.mblend] Joining result into a single raster")) gscript.run_command('r.patch', input=high + ',' + smooth_low_res, output=output) except Exception, ex: gscript.error(_("[r.mblend] ERROR: Failed to create smoothed raster.")) exit() gscript.message(_("[r.mblend] SUCCESS: smoothed raster created.")) if __name__ == '__main__': atexit.register(cleanup) gscript.use_temp_region() main()