1#!/usr/local/bin/python3.8 2 3############################################################################ 4# 5# MODULE: i.pansharpen 6# 7# AUTHOR(S): Overall script by Michael Barton (ASU) 8# Brovey transformation in i.fusion.brovey by Markus Neteler <<neteler at osgeo org>> 9# i.fusion brovey converted to Python by Glynn Clements 10# IHS and PCA transformation added by Michael Barton (ASU) 11# histogram matching algorithm by Michael Barton and Luca Delucchi, Fondazione E. Mach (Italy) 12# Thanks to Markus Metz for help with PCA inversion 13# Thanks to Hamish Bowman for parallel processing algorithm 14# 15# PURPOSE: Sharpening of 3 RGB channels using a high-resolution panchromatic channel 16# 17# COPYRIGHT: (C) 2002-2019 by the GRASS Development Team 18# 19# This program is free software under the GNU General Public 20# License (>=v2). Read the file COPYING that comes with GRASS 21# for details. 22# 23# REFERENCES: 24# Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS and merged MSS/RBV 25# data for analysis of natural vegetation. Proc. of the 14th International 26# Symposium on Remote Sensing of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007. 27# 28# Amarsaikhan, D., & Douglas, T. (2004). Data fusion and multisource image classification. 29# International Journal of Remote Sensing, 25(17), 3529-3539. 30# 31# Behnia, P. (2005). Comparison between four methods for data fusion of ETM+ 32# multispectral and pan images. Geo-spatial Information Science, 8(2), 98-103 33# 34# for LANDSAT 5: see Pohl, C 1996 and others 35# 36############################################################################# 37 38#%Module 39#% description: Image fusion algorithms to sharpen multispectral with high-res panchromatic channels 40#% keyword: imagery 41#% keyword: fusion 42#% keyword: sharpen 43#% keyword: Brovey 44#% keyword: IHS 45#% keyword: HIS 46#% keyword: PCA 47#% overwrite: yes 48#%End 49#%option G_OPT_R_INPUT 50#% key: red 51#% description: Name of raster map to be used for <red> 52#%end 53#%option G_OPT_R_INPUT 54#% key: green 55#% description: Name of raster map to be used for <green> 56#%end 57#%option G_OPT_R_INPUT 58#% key: blue 59#% description: Name of raster map to be used for <blue> 60#%end 61#% option G_OPT_R_INPUT 62#% key: pan 63#% description: Name of raster map to be used for high resolution panchromatic channel 64#%end 65#%option G_OPT_R_BASENAME_OUTPUT 66#%end 67#%option 68#% key: method 69#% description: Method for pan sharpening 70#% options: brovey,ihs,pca 71#% answer: ihs 72#% required: yes 73#%end 74#%option 75#% key: bitdepth 76#% type: integer 77#% description: Bit depth of image (must be in range of 2-30) 78#% options: 2-32 79#% answer: 8 80#% required: yes 81#%end 82#%flag 83#% key: s 84#% description: Serial processing rather than parallel processing 85#%end 86#%flag 87#% key: l 88#% description: Rebalance blue channel for LANDSAT 89#%end 90#%flag 91#% key: r 92#% description: Rescale (stretch) the range of pixel values in each channel to the entire 0-255 8-bit range for processing (see notes) 93#%end 94 95import os 96 97try: 98 import numpy as np 99 hasNumPy = True 100except ImportError: 101 hasNumPy = False 102 103import grass.script as grass 104 105 106def main(): 107 if not hasNumPy: 108 grass.fatal(_("Required dependency NumPy not found. Exiting.")) 109 110 sharpen = options['method'] # sharpening algorithm 111 ms1_orig = options['blue'] # blue channel 112 ms2_orig = options['green'] # green channel 113 ms3_orig = options['red'] # red channel 114 pan_orig = options['pan'] # high res pan channel 115 out = options['output'] # prefix for output RGB maps 116 bits = options['bitdepth'] # bit depth of image channels 117 bladjust = flags['l'] # adjust blue channel 118 sproc = flags['s'] # serial processing 119 rescale = flags['r'] # rescale to spread pixel values to entire 0-255 range 120 121 # Checking bit depth 122 bits = float(bits) 123 if bits < 2 or bits > 30: 124 grass.warning(_("Bit depth is outside acceptable range")) 125 return 126 127 outb = grass.core.find_file('%s_blue' % out) 128 outg = grass.core.find_file('%s_green' % out) 129 outr = grass.core.find_file('%s_red' % out) 130 131 if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and not grass.overwrite(): 132 grass.warning(_('Maps with selected output prefix names already exist.' 133 ' Delete them or use overwrite flag')) 134 return 135 136 pid = str(os.getpid()) 137 138 # convert input image channels to 8 bit for processing 139 ms1 = 'tmp%s_ms1' % pid 140 ms2 = 'tmp%s_ms2' % pid 141 ms3 = 'tmp%s_ms3' % pid 142 pan = 'tmp%s_pan' % pid 143 144 if rescale == False: 145 if bits == 8: 146 grass.message(_("Using 8bit image channels")) 147 if sproc: 148 # serial processing 149 grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1), 150 quiet=True, overwrite=True) 151 grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2), 152 quiet=True, overwrite=True) 153 grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3), 154 quiet=True, overwrite=True) 155 grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan), 156 quiet=True, overwrite=True) 157 else: 158 # parallel processing 159 pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1), 160 quiet=True, overwrite=True) 161 pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2), 162 quiet=True, overwrite=True) 163 pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3), 164 quiet=True, overwrite=True) 165 pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan), 166 quiet=True, overwrite=True) 167 168 pb.wait() 169 pg.wait() 170 pr.wait() 171 pp.wait() 172 173 else: 174 grass.message(_("Converting image chanels to 8bit for processing")) 175 maxval = pow(2, bits) - 1 176 if sproc: 177 # serial processing 178 grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval, 179 output=ms1, to='0,255', quiet=True, overwrite=True) 180 grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval, 181 output=ms2, to='0,255', quiet=True, overwrite=True) 182 grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval, 183 output=ms3, to='0,255', quiet=True, overwrite=True) 184 grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval, 185 output=pan, to='0,255', quiet=True, overwrite=True) 186 187 else: 188 # parallel processing 189 pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval, 190 output=ms1, to='0,255', quiet=True, overwrite=True) 191 pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval, 192 output=ms2, to='0,255', quiet=True, overwrite=True) 193 pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval, 194 output=ms3, to='0,255', quiet=True, overwrite=True) 195 pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval, 196 output=pan, to='0,255', quiet=True, overwrite=True) 197 198 pb.wait() 199 pg.wait() 200 pr.wait() 201 pp.wait() 202 203 else: 204 grass.message(_("Rescaling image chanels to 8bit for processing")) 205 206 min_ms1 = int(grass.raster_info(ms1_orig)['min']) 207 max_ms1 = int(grass.raster_info(ms1_orig)['max']) 208 min_ms2 = int(grass.raster_info(ms2_orig)['min']) 209 max_ms2 = int(grass.raster_info(ms2_orig)['max']) 210 min_ms3 = int(grass.raster_info(ms3_orig)['min']) 211 max_ms3 = int(grass.raster_info(ms3_orig)['max']) 212 min_pan = int(grass.raster_info(pan_orig)['min']) 213 max_pan = int(grass.raster_info(pan_orig)['max']) 214 215 maxval = pow(2, bits) - 1 216 if sproc: 217 # serial processing 218 grass.run_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1), 219 output=ms1, to='0,255', quiet=True, overwrite=True) 220 grass.run_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2), 221 output=ms2, to='0,255', quiet=True, overwrite=True) 222 grass.run_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3), 223 output=ms3, to='0,255', quiet=True, overwrite=True) 224 grass.run_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan), 225 output=pan, to='0,255', quiet=True, overwrite=True) 226 227 else: 228 # parallel processing 229 pb = grass.start_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1), 230 output=ms1, to='0,255', quiet=True, overwrite=True) 231 pg = grass.start_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2), 232 output=ms2, to='0,255', quiet=True, overwrite=True) 233 pr = grass.start_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3), 234 output=ms3, to='0,255', quiet=True, overwrite=True) 235 pp = grass.start_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan), 236 output=pan, to='0,255', quiet=True, overwrite=True) 237 238 pb.wait() 239 pg.wait() 240 pr.wait() 241 pp.wait() 242 243 244 # get PAN resolution: 245 kv = grass.raster_info(map=pan) 246 nsres = kv['nsres'] 247 ewres = kv['ewres'] 248 panres = (nsres + ewres) / 2 249 250 # clone current region 251 grass.use_temp_region() 252 grass.run_command('g.region', res=panres, align=pan) 253 254 # Select sharpening method 255 grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres)) 256 if sharpen == "brovey": 257 brovey(pan, ms1, ms2, ms3, out, pid, sproc) 258 elif sharpen == "ihs": 259 ihs(pan, ms1, ms2, ms3, out, pid, sproc) 260 elif sharpen == "pca": 261 pca(pan, ms1, ms2, ms3, out, pid, sproc) 262 # Could add other sharpening algorithms here, e.g. wavelet transformation 263 264 grass.message(_("Assigning grey equalized color tables to output images...")) 265 266 # equalized grey scales give best contrast 267 grass.message(_("setting pan-sharpened channels to equalized grey scale")) 268 for ch in ['red', 'green', 'blue']: 269 grass.run_command('r.colors', quiet=True, map="%s_%s" % (out, ch), 270 flags="e", color='grey') 271 272 # Landsat too blue-ish because panchromatic band less sensitive to blue 273 # light, so output blue channed can be modified 274 if bladjust: 275 grass.message(_("Adjusting blue channel color table...")) 276 blue_colors = ['0 0 0 0\n5% 0 0 0\n67% 255 255 255\n100% 255 255 255'] 277 # these previous colors are way too blue for landsat 278 # blue_colors = ['0 0 0 0\n10% 0 0 0\n20% 200 200 200\n40% 230 230 230\n67% 255 255 255\n100% 255 255 255'] 279 bc = grass.feed_command('r.colors', quiet = True, map = "%s_blue" % out, rules = "-") 280 bc.stdin.write(grass.encode('\n'.join(blue_colors))) 281 bc.stdin.close() 282 283 # output notice 284 grass.verbose(_("The following pan-sharpened output maps have been generated:")) 285 for ch in ['red', 'green', 'blue']: 286 grass.verbose(_("%s_%s") % (out, ch)) 287 288 grass.verbose(_("To visualize output, run: g.region -p raster=%s_red" % out)) 289 grass.verbose(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out))) 290 grass.verbose(_("If desired, combine channels into a single RGB map with 'r.composite'.")) 291 grass.verbose(_("Channel colors can be rebalanced using i.colors.enhance.")) 292 293 # write cmd history: 294 for ch in ['red', 'green', 'blue']: 295 grass.raster_history("%s_%s" % (out, ch)) 296 297 # create a group with the three outputs 298 #grass.run_command('i.group', group=out, 299 # input="{n}_red,{n}_blue,{n}_green".format(n=out)) 300 301 # Cleanup 302 grass.message(_("cleaning up temp files")) 303 try: 304 grass.run_command('g.remove', flags="f", type="raster", 305 pattern="tmp%s*" % pid, quiet=True) 306 except: 307 "" 308 309def brovey(pan, ms1, ms2, ms3, out, pid, sproc): 310 grass.verbose(_("Using Brovey algorithm")) 311 312 # pan/intensity histogram matching using linear regression 313 grass.message(_("Pan channel/intensity histogram matching using linear regression")) 314 outname = 'tmp%s_pan1' % pid 315 panmatch1 = matchhist(pan, ms1, outname) 316 317 outname = 'tmp%s_pan2' % pid 318 panmatch2 = matchhist(pan, ms2, outname) 319 320 outname = 'tmp%s_pan3' % pid 321 panmatch3 = matchhist(pan, ms3, outname) 322 323 outr = '%s_red' % out 324 outg = '%s_green' % out 325 outb = '%s_blue' % out 326 327 # calculate brovey transformation 328 grass.message(_("Calculating Brovey transformation...")) 329 330 if sproc: 331 # serial processing 332 e = '''eval(k = "$ms1" + "$ms2" + "$ms3") 333 "$outr" = 1 * round("$ms3" * "$panmatch3" / k) 334 "$outg" = 1 * round("$ms2" * "$panmatch2" / k) 335 "$outb" = 1 * round("$ms1" * "$panmatch1" / k)''' 336 grass.mapcalc(e, outr=outr, outg=outg, outb=outb, 337 panmatch1=panmatch1, panmatch2=panmatch2, 338 panmatch3=panmatch3, ms1=ms1, ms2=ms2, ms3=ms3, 339 overwrite=True) 340 else: 341 # parallel processing 342 pb = grass.mapcalc_start('%s_blue = 1 * round((%s * %s) / (%s + %s + %s))' % 343 (out, ms1, panmatch1, ms1, ms2, ms3), 344 overwrite=True) 345 pg = grass.mapcalc_start('%s_green = 1 * round((%s * %s) / (%s + %s + %s))' % 346 (out, ms2, panmatch2, ms1, ms2, ms3), 347 overwrite=True) 348 pr = grass.mapcalc_start('%s_red = 1 * round((%s * %s) / (%s + %s + %s))' % 349 (out, ms3, panmatch3, ms1, ms2, ms3), 350 overwrite=True) 351 352 pb.wait(), pg.wait(), pr.wait() 353 try: 354 pb.terminate(), pg.terminate(), pr.terminate() 355 except: 356 "" 357 358 # Cleanup 359 try: 360 grass.run_command('g.remove', flags='f', quiet=True, type='raster', 361 name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3)) 362 except: 363 "" 364 365def ihs(pan, ms1, ms2, ms3, out, pid, sproc): 366 grass.verbose(_("Using IHS<->RGB algorithm")) 367 # transform RGB channels into IHS color space 368 grass.message(_("Transforming to IHS color space...")) 369 grass.run_command('i.rgb.his', overwrite=True, 370 red=ms3, 371 green=ms2, 372 blue=ms1, 373 hue="tmp%s_hue" % pid, 374 intensity="tmp%s_int" % pid, 375 saturation="tmp%s_sat" % pid) 376 377 # pan/intensity histogram matching using linear regression 378 target = "tmp%s_int" % pid 379 outname = "tmp%s_pan_int" % pid 380 panmatch = matchhist(pan, target, outname) 381 382 # substitute pan for intensity channel and transform back to RGB color space 383 grass.message(_("Transforming back to RGB color space and sharpening...")) 384 grass.run_command('i.his.rgb', overwrite=True, 385 hue="tmp%s_hue" % pid, 386 intensity="%s" % panmatch, 387 saturation="tmp%s_sat" % pid, 388 red="%s_red" % out, 389 green="%s_green" % out, 390 blue="%s_blue" % out) 391 392 # Cleanup 393 try: 394 grass.run_command('g.remove', flags='f', quiet=True, type='raster', 395 name=panmatch) 396 except: 397 "" 398 399def pca(pan, ms1, ms2, ms3, out, pid, sproc): 400 401 grass.verbose(_("Using PCA/inverse PCA algorithm")) 402 grass.message(_("Creating PCA images and calculating eigenvectors...")) 403 404 # initial PCA with RGB channels 405 pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0', 406 input='%s,%s,%s' % (ms1, ms2, ms3), 407 output='tmp%s.pca' % pid) 408 if len(pca_out) < 1: 409 grass.fatal(_("Input has no data. Check region settings.")) 410 411 b1evect = [] 412 b2evect = [] 413 b3evect = [] 414 for l in pca_out.replace('(', ',').replace(')', ',').splitlines(): 415 b1evect.append(float(l.split(',')[1])) 416 b2evect.append(float(l.split(',')[2])) 417 b3evect.append(float(l.split(',')[3])) 418 419 # inverse PCA with hi res pan channel substituted for principal component 1 420 pca1 = 'tmp%s.pca.1' % pid 421 pca2 = 'tmp%s.pca.2' % pid 422 pca3 = 'tmp%s.pca.3' % pid 423 b1evect1 = b1evect[0] 424 b1evect2 = b1evect[1] 425 b1evect3 = b1evect[2] 426 b2evect1 = b2evect[0] 427 b2evect2 = b2evect[1] 428 b2evect3 = b2evect[2] 429 b3evect1 = b3evect[0] 430 b3evect2 = b3evect[1] 431 b3evect3 = b3evect[2] 432 433 # Histogram matching 434 outname = 'tmp%s_pan1' % pid 435 panmatch1 = matchhist(pan, ms1, outname) 436 437 outname = 'tmp%s_pan2' % pid 438 panmatch2 = matchhist(pan, ms2, outname) 439 440 outname = 'tmp%s_pan3' % pid 441 panmatch3 = matchhist(pan, ms3, outname) 442 443 grass.message(_("Performing inverse PCA ...")) 444 445 # Get mean value of each channel 446 stats1 = grass.parse_command("r.univar", map=ms1, flags='g', 447 parse=(grass.parse_key_val, 448 {'sep': '='})) 449 stats2 = grass.parse_command("r.univar", map=ms2, flags='g', 450 parse=(grass.parse_key_val, 451 {'sep': '='})) 452 stats3 = grass.parse_command("r.univar", map=ms3, flags='g', 453 parse=(grass.parse_key_val, 454 {'sep': '='})) 455 456 b1mean = float(stats1['mean']) 457 b2mean = float(stats2['mean']) 458 b3mean = float(stats3['mean']) 459 460 if sproc: 461 # serial processing 462 outr = '%s_red' % out 463 outg = '%s_green' % out 464 outb = '%s_blue' % out 465 466 cmd1 = "$outb = 1 * round(($panmatch1 * $b1evect1) + ($pca2 * $b1evect2) + ($pca3 * $b1evect3) + $b1mean)" 467 cmd2 = "$outg = 1 * round(($panmatch2 * $b2evect1) + ($pca2 * $b2evect2) + ($pca3 * $b2evect3) + $b2mean)" 468 cmd3 = "$outr = 1 * round(($panmatch3 * $b3evect1) + ($pca2 * $b3evect2) + ($pca3 * $b3evect3) + $b3mean)" 469 470 cmd = '\n'.join([cmd1, cmd2, cmd3]) 471 472 grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr, 473 panmatch1=panmatch1, 474 panmatch2=panmatch2, 475 panmatch3=panmatch3, 476 pca2=pca2, 477 pca3=pca3, 478 b1evect1=b1evect1, 479 b2evect1=b2evect1, 480 b3evect1=b3evect1, 481 b1evect2=b1evect2, 482 b2evect2=b2evect2, 483 b3evect2=b3evect2, 484 b1evect3=b1evect3, 485 b2evect3=b2evect3, 486 b3evect3=b3evect3, 487 b1mean=b1mean, 488 b2mean=b2mean, 489 b3mean=b3mean, 490 overwrite=True) 491 else: 492 # parallel processing 493 pb = grass.mapcalc_start('%s_blue = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)' 494 % (out, 495 panmatch1, 496 b1evect1, 497 pca2, 498 b1evect2, 499 pca3, 500 b1evect3, 501 b1mean), 502 overwrite=True) 503 504 pg = grass.mapcalc_start('%s_green = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)' 505 % (out, 506 panmatch2, 507 b2evect1, 508 pca2, 509 b2evect2, 510 pca3, 511 b2evect3, 512 b2mean), 513 overwrite=True) 514 515 pr = grass.mapcalc_start('%s_red = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)' 516 % (out, 517 panmatch3, 518 b3evect1, 519 pca2, 520 b3evect2, 521 pca3, 522 b3evect3, 523 b3mean), 524 overwrite=True) 525 526 pb.wait(), pg.wait(), pr.wait() 527 try: 528 pb.terminate(), pg.terminate(), pr.terminate() 529 except: 530 "" 531 532 # Cleanup 533 grass.run_command('g.remove', flags='f', quiet=True, type='raster', 534 name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3)) 535 536def matchhist(original, target, matched): 537 # pan/intensity histogram matching using numpy arrays 538 grass.message(_("Histogram matching...")) 539 540 # input images 541 original = original.split('@')[0] 542 target = target.split('@')[0] 543 images = [original, target] 544 545 # create a dictionary to hold arrays for each image 546 arrays = {} 547 548 for img in images: 549 # calculate number of cells for each grey value for for each image 550 stats_out = grass.pipe_command('r.stats', flags='cin', input=img, 551 sep=':') 552 stats = grass.decode(stats_out.communicate()[0]).split('\n')[:-1] 553 stats_dict = dict(s.split(':', 1) for s in stats) 554 total_cells = 0 # total non-null cells 555 for j in stats_dict: 556 stats_dict[j] = int(stats_dict[j]) 557 if j != '*': 558 total_cells += stats_dict[j] 559 560 if total_cells < 1: 561 grass.fatal(_("Input has no data. Check region settings.")) 562 563 # Make a 2x256 structured array for each image with a 564 # cumulative distribution function (CDF) for each grey value. 565 # Grey value is the integer (i4) and cdf is float (f4). 566 567 arrays[img] = np.zeros((256, ), dtype=('i4,f4')) 568 cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values 569 570 for n in range(0, 256): 571 if str(n) in stats_dict: 572 num_cells = stats_dict[str(n)] 573 else: 574 num_cells = 0 575 576 cum_cells += num_cells 577 578 # cdf is the the number of cells at or below a given grey value 579 # divided by the total number of cells 580 cdf = float(cum_cells) / float(total_cells) 581 582 # insert values into array 583 arrays[img][n] = (n, cdf) 584 585 # open file for reclass rules 586 outfile = open(grass.tempfile(), 'w') 587 588 for i in arrays[original]: 589 # for each grey value and corresponding cdf value in original, find the 590 # cdf value in target that is closest to the target cdf value 591 difference_list = [] 592 for j in arrays[target]: 593 # make a list of the difference between each original cdf value and 594 # the target cdf value 595 difference_list.append(abs(i[1] - j[1])) 596 597 # get the smallest difference in the list 598 min_difference = min(difference_list) 599 600 for j in arrays[target]: 601 # find the grey value in target that corresponds to the cdf 602 # closest to the original cdf 603 if j[1] <= i[1] + min_difference and j[1] >= i[1] - min_difference: 604 # build a reclass rules file from the original grey value and 605 # corresponding grey value from target 606 out_line = "%d = %d\n" % (i[0], j[0]) 607 outfile.write(out_line) 608 break 609 610 outfile.close() 611 612 # create reclass of target from reclass rules file 613 result = grass.core.find_file(matched, element='cell') 614 if result['fullname']: 615 grass.run_command('g.remove', flags='f', quiet=True, type='raster', 616 name=matched) 617 grass.run_command('r.reclass', input=original, out=matched, 618 rules=outfile.name) 619 else: 620 grass.run_command('r.reclass', input=original, out=matched, 621 rules=outfile.name) 622 623 # Cleanup 624 # remove the rules file 625 grass.try_remove(outfile.name) 626 627 # return reclass of target with histogram that matches original 628 return matched 629 630if __name__ == "__main__": 631 options, flags = grass.parser() 632 main() 633