| 1 | #!/usr/bin/env python
|
|---|
| 2 |
|
|---|
| 3 | ############################################################################
|
|---|
| 4 | #
|
|---|
| 5 | # MODULE: i.panmethod
|
|---|
| 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-2012 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 | #%flag
|
|---|
| 75 | #% key: s
|
|---|
| 76 | #% description: Serial processing rather than parallel processing
|
|---|
| 77 | #%end
|
|---|
| 78 | #%flag
|
|---|
| 79 | #% key: l
|
|---|
| 80 | #% description: Rebalance blue channel for LANDSAT
|
|---|
| 81 | #%end
|
|---|
| 82 |
|
|---|
| 83 | import os
|
|---|
| 84 |
|
|---|
| 85 | try:
|
|---|
| 86 | import numpy as np
|
|---|
| 87 | hasNumPy = True
|
|---|
| 88 | except ImportError:
|
|---|
| 89 | hasNumPy = False
|
|---|
| 90 |
|
|---|
| 91 | import grass.script as grass
|
|---|
| 92 |
|
|---|
| 93 |
|
|---|
| 94 | def main():
|
|---|
| 95 | if not hasNumPy:
|
|---|
| 96 | grass.fatal(_("Required dependency NumPy not found. Exiting."))
|
|---|
| 97 |
|
|---|
| 98 | sharpen = options['method'] # sharpening algorithm
|
|---|
| 99 | ms1 = options['blue'] # blue channel
|
|---|
| 100 | ms2 = options['green'] # green channel
|
|---|
| 101 | ms3 = options['red'] # red channel
|
|---|
| 102 | pan = options['pan'] # high res pan channel
|
|---|
| 103 | out = options['output'] # prefix for output RGB maps
|
|---|
| 104 | bladjust = flags['l'] # adjust blue channel
|
|---|
| 105 | sproc = flags['s'] # serial processing
|
|---|
| 106 |
|
|---|
| 107 | outb = grass.core.find_file('%s_blue' % out)
|
|---|
| 108 | outg = grass.core.find_file('%s_green' % out)
|
|---|
| 109 | outr = grass.core.find_file('%s_red' % out)
|
|---|
| 110 |
|
|---|
| 111 | if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and not grass.overwrite():
|
|---|
| 112 | grass.warning(_('Maps with selected output prefix names already exist.'
|
|---|
| 113 | ' Delete them or use overwrite flag'))
|
|---|
| 114 | return
|
|---|
| 115 |
|
|---|
| 116 | pid = str(os.getpid())
|
|---|
| 117 |
|
|---|
| 118 | # get PAN resolution:
|
|---|
| 119 | kv = grass.raster_info(map=pan)
|
|---|
| 120 | nsres = kv['nsres']
|
|---|
| 121 | ewres = kv['ewres']
|
|---|
| 122 | panres = (nsres + ewres) / 2
|
|---|
| 123 |
|
|---|
| 124 | # clone current region
|
|---|
| 125 | grass.use_temp_region()
|
|---|
| 126 |
|
|---|
| 127 | grass.run_command('g.region', res=panres, align=pan)
|
|---|
| 128 |
|
|---|
| 129 | grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
|
|---|
| 130 |
|
|---|
| 131 | if sharpen == "brovey":
|
|---|
| 132 | grass.verbose(_("Using Brovey algorithm"))
|
|---|
| 133 |
|
|---|
| 134 | # pan/intensity histogram matching using linear regression
|
|---|
| 135 | outname = 'tmp%s_pan1' % pid
|
|---|
| 136 | panmatch1 = matchhist(pan, ms1, outname)
|
|---|
| 137 |
|
|---|
| 138 | outname = 'tmp%s_pan2' % pid
|
|---|
| 139 | panmatch2 = matchhist(pan, ms2, outname)
|
|---|
| 140 |
|
|---|
| 141 | outname = 'tmp%s_pan3' % pid
|
|---|
| 142 | panmatch3 = matchhist(pan, ms3, outname)
|
|---|
| 143 |
|
|---|
| 144 | outr = '%s_red' % out
|
|---|
| 145 | outg = '%s_green' % out
|
|---|
| 146 | outb = '%s_blue' % out
|
|---|
| 147 |
|
|---|
| 148 | # calculate brovey transformation
|
|---|
| 149 | grass.message(_("Calculating Brovey transformation..."))
|
|---|
| 150 |
|
|---|
| 151 | if sproc:
|
|---|
| 152 | # serial processing
|
|---|
| 153 | e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
|
|---|
| 154 | "$outr" = 1.0 * "$ms3" * "$panmatch3" / k
|
|---|
| 155 | "$outg" = 1.0 * "$ms2" * "$panmatch2" / k
|
|---|
| 156 | "$outb" = 1.0 * "$ms1" * "$panmatch1" / k'''
|
|---|
| 157 | grass.mapcalc(e, outr=outr, outg=outg, outb=outb,
|
|---|
| 158 | panmatch1=panmatch1, panmatch2=panmatch2,
|
|---|
| 159 | panmatch3=panmatch3, ms1=ms1, ms2=ms2, ms3=ms3,
|
|---|
| 160 | overwrite=True)
|
|---|
| 161 | else:
|
|---|
| 162 | # parallel processing
|
|---|
| 163 | pb = grass.mapcalc_start('%s_blue = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 164 | (out, ms1, panmatch1, ms1, ms2, ms3),
|
|---|
| 165 | overwrite=True)
|
|---|
| 166 | pg = grass.mapcalc_start('%s_green = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 167 | (out, ms2, panmatch2, ms1, ms2, ms3),
|
|---|
| 168 | overwrite=True)
|
|---|
| 169 | pr = grass.mapcalc_start('%s_red = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 170 | (out, ms3, panmatch3, ms1, ms2, ms3),
|
|---|
| 171 | overwrite=True)
|
|---|
| 172 |
|
|---|
| 173 | pb.wait()
|
|---|
| 174 | pg.wait()
|
|---|
| 175 | pr.wait()
|
|---|
| 176 |
|
|---|
| 177 | # Cleanup
|
|---|
| 178 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 179 | name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|---|
| 180 |
|
|---|
| 181 | elif sharpen == "ihs":
|
|---|
| 182 | grass.verbose(_("Using IHS<->RGB algorithm"))
|
|---|
| 183 | # transform RGB channels into IHS color space
|
|---|
| 184 | grass.message(_("Transforming to IHS color space..."))
|
|---|
| 185 | grass.run_command('i.rgb.his', overwrite=True,
|
|---|
| 186 | red=ms3,
|
|---|
| 187 | green=ms2,
|
|---|
| 188 | blue=ms1,
|
|---|
| 189 | hue="tmp%s_hue" % pid,
|
|---|
| 190 | intensity="tmp%s_int" % pid,
|
|---|
| 191 | saturation="tmp%s_sat" % pid)
|
|---|
| 192 |
|
|---|
| 193 | # pan/intensity histogram matching using linear regression
|
|---|
| 194 | target = "tmp%s_int" % pid
|
|---|
| 195 | outname = "tmp%s_pan_int" % pid
|
|---|
| 196 | panmatch = matchhist(pan, target, outname)
|
|---|
| 197 |
|
|---|
| 198 | # substitute pan for intensity channel and transform back to RGB color space
|
|---|
| 199 | grass.message(_("Transforming back to RGB color space and sharpening..."))
|
|---|
| 200 | grass.run_command('i.his.rgb', overwrite=True,
|
|---|
| 201 | hue="tmp%s_hue" % pid,
|
|---|
| 202 | intensity="%s" % panmatch,
|
|---|
| 203 | saturation="tmp%s_sat" % pid,
|
|---|
| 204 | red="%s_red" % out,
|
|---|
| 205 | green="%s_green" % out,
|
|---|
| 206 | blue="%s_blue" % out)
|
|---|
| 207 |
|
|---|
| 208 | # Cleanup
|
|---|
| 209 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 210 | name=panmatch)
|
|---|
| 211 |
|
|---|
| 212 | elif sharpen == "pca":
|
|---|
| 213 | grass.verbose(_("Using PCA/inverse PCA algorithm"))
|
|---|
| 214 | grass.message(_("Creating PCA images and calculating eigenvectors..."))
|
|---|
| 215 |
|
|---|
| 216 | # initial PCA with RGB channels
|
|---|
| 217 | pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0',
|
|---|
| 218 | input='%s,%s,%s' % (ms1, ms2, ms3),
|
|---|
| 219 | output='tmp%s.pca' % pid)
|
|---|
| 220 | if len(pca_out) < 1:
|
|---|
| 221 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 222 |
|
|---|
| 223 | b1evect = []
|
|---|
| 224 | b2evect = []
|
|---|
| 225 | b3evect = []
|
|---|
| 226 | for l in pca_out.replace('(', ',').replace(')', ',').splitlines():
|
|---|
| 227 | b1evect.append(float(l.split(',')[1]))
|
|---|
| 228 | b2evect.append(float(l.split(',')[2]))
|
|---|
| 229 | b3evect.append(float(l.split(',')[3]))
|
|---|
| 230 |
|
|---|
| 231 | # inverse PCA with hi res pan channel substituted for principal component 1
|
|---|
| 232 | pca1 = 'tmp%s.pca.1' % pid
|
|---|
| 233 | pca2 = 'tmp%s.pca.2' % pid
|
|---|
| 234 | pca3 = 'tmp%s.pca.3' % pid
|
|---|
| 235 | b1evect1 = b1evect[0]
|
|---|
| 236 | b1evect2 = b1evect[1]
|
|---|
| 237 | b1evect3 = b1evect[2]
|
|---|
| 238 | b2evect1 = b2evect[0]
|
|---|
| 239 | b2evect2 = b2evect[1]
|
|---|
| 240 | b2evect3 = b2evect[2]
|
|---|
| 241 | b3evect1 = b3evect[0]
|
|---|
| 242 | b3evect2 = b3evect[1]
|
|---|
| 243 | b3evect3 = b3evect[2]
|
|---|
| 244 |
|
|---|
| 245 | outname = 'tmp%s_pan' % pid
|
|---|
| 246 | panmatch = matchhist(pan, ms1, outname)
|
|---|
| 247 |
|
|---|
| 248 | grass.message(_("Performing inverse PCA ..."))
|
|---|
| 249 |
|
|---|
| 250 | stats1 = grass.parse_command("r.univar", map=ms1, flags='g',
|
|---|
| 251 | parse=(grass.parse_key_val,
|
|---|
| 252 | {'sep': '='}))
|
|---|
| 253 | stats2 = grass.parse_command("r.univar", map=ms2, flags='g',
|
|---|
| 254 | parse=(grass.parse_key_val,
|
|---|
| 255 | {'sep': '='}))
|
|---|
| 256 | stats3 = grass.parse_command("r.univar", map=ms3, flags='g',
|
|---|
| 257 | parse=(grass.parse_key_val,
|
|---|
| 258 | {'sep': '='}))
|
|---|
| 259 |
|
|---|
| 260 | b1mean = float(stats1['mean'])
|
|---|
| 261 | b2mean = float(stats2['mean'])
|
|---|
| 262 | b3mean = float(stats3['mean'])
|
|---|
| 263 |
|
|---|
| 264 | if sproc:
|
|---|
| 265 | # serial processing
|
|---|
| 266 | e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
|
|---|
| 267 | "$outr" = 1.0 * "$ms3" * "$panmatch3" / k
|
|---|
| 268 | "$outg" = 1.0 * "$ms2" * "$panmatch2" / k
|
|---|
| 269 | "$outb" = 1.0* "$ms1" * "$panmatch1" / k'''
|
|---|
| 270 |
|
|---|
| 271 | outr = '%s_red' % out
|
|---|
| 272 | outg = '%s_green' % out
|
|---|
| 273 | outb = '%s_blue' % out
|
|---|
| 274 |
|
|---|
| 275 | cmd1 = "$outb = (1.0 * $panmatch * $b1evect1) + ($pca2 * $b2evect1) + ($pca3 * $b3evect1) + $b1mean"
|
|---|
| 276 | cmd2 = "$outg = (1.0 * $panmatch * $b1evect2) + ($pca2 * $b2evect1) + ($pca3 * $b3evect2) + $b2mean"
|
|---|
| 277 | cmd3 = "$outr = (1.0 * $panmatch * $b1evect3) + ($pca2 * $b2evect3) + ($pca3 * $b3evect3) + $b3mean"
|
|---|
| 278 |
|
|---|
| 279 | cmd = '\n'.join([cmd1, cmd2, cmd3])
|
|---|
| 280 |
|
|---|
| 281 | grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr,
|
|---|
| 282 | panmatch=panmatch, pca2=pca2, pca3=pca3,
|
|---|
| 283 | b1evect1=b1evect1, b2evect1=b2evect1, b3evect1=b3evect1,
|
|---|
| 284 | b1evect2=b1evect2, b2evect2=b2evect2, b3evect2=b3evect2,
|
|---|
| 285 | b1evect3=b1evect3, b2evect3=b2evect3, b3evect3=b3evect3,
|
|---|
| 286 | b1mean=b1mean, b2mean=b2mean, b3mean=b3mean,
|
|---|
| 287 | overwrite=True)
|
|---|
| 288 | else:
|
|---|
| 289 | # parallel processing
|
|---|
| 290 | pb = grass.mapcalc_start('%s_blue = (%s * %f) + (%s * %f) + (%s * %f) + %f'
|
|---|
| 291 | % (out, panmatch, b1evect1, pca2,
|
|---|
| 292 | b2evect1, pca3, b3evect1, b1mean),
|
|---|
| 293 | overwrite=True)
|
|---|
| 294 |
|
|---|
| 295 | pg = grass.mapcalc_start('%s_green = (%s * %f) + (%s * %f) + (%s * %f) + %f'
|
|---|
| 296 | % (out, panmatch, b1evect2, pca2,
|
|---|
| 297 | b2evect2, pca3, b3evect2, b2mean),
|
|---|
| 298 | overwrite=True)
|
|---|
| 299 |
|
|---|
| 300 | pr = grass.mapcalc_start('%s_red = (%s * %f) + (%s * %f) + (%s * ''%f) + %f'
|
|---|
| 301 | % (out, panmatch, b1evect3, pca2,
|
|---|
| 302 | b2evect3, pca3, b3evect3, b3mean),
|
|---|
| 303 | overwrite=True)
|
|---|
| 304 |
|
|---|
| 305 | pr.wait()
|
|---|
| 306 | pg.wait()
|
|---|
| 307 | pb.wait()
|
|---|
| 308 |
|
|---|
| 309 | # Cleanup
|
|---|
| 310 | grass.run_command('g.remove', flags='f', quiet=True, type="raster",
|
|---|
| 311 | pattern='tmp%s*,%s' % (pid, panmatch))
|
|---|
| 312 |
|
|---|
| 313 | # Could add other sharpening algorithms here, e.g. wavelet transformation
|
|---|
| 314 |
|
|---|
| 315 | grass.message(_("Assigning grey equalized color tables to output images..."))
|
|---|
| 316 | # equalized grey scales give best contrast
|
|---|
| 317 | for ch in ['red', 'green', 'blue']:
|
|---|
| 318 | grass.run_command('r.colors', quiet=True, map="%s_%s" % (out, ch),
|
|---|
| 319 | flags="e", color='grey')
|
|---|
| 320 |
|
|---|
| 321 | # Landsat too blue-ish because panchromatic band less sensitive to blue
|
|---|
| 322 | # light, so output blue channed can be modified
|
|---|
| 323 | if bladjust:
|
|---|
| 324 | grass.message(_("Adjusting blue channel color table..."))
|
|---|
| 325 | rules = grass.tempfile()
|
|---|
| 326 | colors = open(rules, 'w')
|
|---|
| 327 | colors.write('5 0 0 0\n20 200 200 200\n40 230 230 230\n67 255 255 255 \n')
|
|---|
| 328 | colors.close()
|
|---|
| 329 |
|
|---|
| 330 | grass.run_command('r.colors', map="%s_blue" % out, rules=rules)
|
|---|
| 331 | os.remove(rules)
|
|---|
| 332 |
|
|---|
| 333 | # output notice
|
|---|
| 334 | grass.verbose(_("The following pan-sharpened output maps have been generated:"))
|
|---|
| 335 | for ch in ['red', 'green', 'blue']:
|
|---|
| 336 | grass.verbose(_("%s_%s") % (out, ch))
|
|---|
| 337 |
|
|---|
| 338 | grass.verbose(_("To visualize output, run: g.region -p raster=%s_red" % out))
|
|---|
| 339 | grass.verbose(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out)))
|
|---|
| 340 | grass.verbose(_("If desired, combine channels into a single RGB map with 'r.composite'."))
|
|---|
| 341 | grass.verbose(_("Channel colors can be rebalanced using i.colors.enhance."))
|
|---|
| 342 |
|
|---|
| 343 | # write cmd history:
|
|---|
| 344 | for ch in ['red', 'green', 'blue']:
|
|---|
| 345 | grass.raster_history("%s_%s" % (out, ch))
|
|---|
| 346 |
|
|---|
| 347 | # create a group with the three output
|
|---|
| 348 | grass.run_command('i.group', group=out,
|
|---|
| 349 | input="{n}_red,{n}_blue,{n}_green".format(n=out))
|
|---|
| 350 |
|
|---|
| 351 | # Cleanup
|
|---|
| 352 | grass.run_command('g.remove', flags="f", type="raster",
|
|---|
| 353 | pattern="tmp%s*" % pid, quiet=True)
|
|---|
| 354 |
|
|---|
| 355 |
|
|---|
| 356 | def matchhist(original, target, matched):
|
|---|
| 357 | # pan/intensity histogram matching using numpy arrays
|
|---|
| 358 | grass.message(_("Histogram matching..."))
|
|---|
| 359 |
|
|---|
| 360 | # input images
|
|---|
| 361 | original = original.split('@')[0]
|
|---|
| 362 | target = target.split('@')[0]
|
|---|
| 363 | images = [original, target]
|
|---|
| 364 |
|
|---|
| 365 | # create a dictionary to hold arrays for each image
|
|---|
| 366 | arrays = {}
|
|---|
| 367 |
|
|---|
| 368 | for i in images:
|
|---|
| 369 | # calculate number of cells for each grey value for for each image
|
|---|
| 370 | stats_out = grass.pipe_command('r.stats', flags='cin', input=i,
|
|---|
| 371 | sep=':')
|
|---|
| 372 | stats = stats_out.communicate()[0].split('\n')[:-1]
|
|---|
| 373 | stats_dict = dict(s.split(':', 1) for s in stats)
|
|---|
| 374 | total_cells = 0 # total non-null cells
|
|---|
| 375 | for j in stats_dict:
|
|---|
| 376 | stats_dict[j] = int(stats_dict[j])
|
|---|
| 377 | if j != '*':
|
|---|
| 378 | total_cells += stats_dict[j]
|
|---|
| 379 |
|
|---|
| 380 | if total_cells < 1:
|
|---|
| 381 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 382 |
|
|---|
| 383 | # Make a 2x256 structured array for each image with a
|
|---|
| 384 | # cumulative distribution function (CDF) for each grey value.
|
|---|
| 385 | # Grey value is the integer (i4) and cdf is float (f4).
|
|---|
| 386 |
|
|---|
| 387 | arrays[i] = np.zeros((256, ), dtype=('i4,f4'))
|
|---|
| 388 | cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values
|
|---|
| 389 |
|
|---|
| 390 | for n in range(0, 256):
|
|---|
| 391 | if str(n) in stats_dict:
|
|---|
| 392 | num_cells = stats_dict[str(n)]
|
|---|
| 393 | else:
|
|---|
| 394 | num_cells = 0
|
|---|
| 395 |
|
|---|
| 396 | cum_cells += num_cells
|
|---|
| 397 |
|
|---|
| 398 | # cdf is the the number of cells at or below a given grey value
|
|---|
| 399 | # divided by the total number of cells
|
|---|
| 400 | cdf = float(cum_cells) / float(total_cells)
|
|---|
| 401 |
|
|---|
| 402 | # insert values into array
|
|---|
| 403 | arrays[i][n] = (n, cdf)
|
|---|
| 404 |
|
|---|
| 405 | # open file for reclass rules
|
|---|
| 406 | outfile = open(grass.tempfile(), 'w')
|
|---|
| 407 |
|
|---|
| 408 | for i in arrays[original]:
|
|---|
| 409 | # for each grey value and corresponding cdf value in original, find the
|
|---|
| 410 | # cdf value in target that is closest to the target cdf value
|
|---|
| 411 | difference_list = []
|
|---|
| 412 | for j in arrays[target]:
|
|---|
| 413 | # make a list of the difference between each original cdf value and
|
|---|
| 414 | # the target cdf value
|
|---|
| 415 | difference_list.append(abs(i[1] - j[1]))
|
|---|
| 416 |
|
|---|
| 417 | # get the smallest difference in the list
|
|---|
| 418 | min_difference = min(difference_list)
|
|---|
| 419 |
|
|---|
| 420 | for j in arrays[target]:
|
|---|
| 421 | # find the grey value in target that corresponds to the cdf
|
|---|
| 422 | # closest to the original cdf
|
|---|
| 423 | if j[1] == i[1] + min_difference or j[1] == i[1] - min_difference:
|
|---|
| 424 | # build a reclass rules file from the original grey value and
|
|---|
| 425 | # corresponding grey value from target
|
|---|
| 426 | out_line = "%d = %d\n" % (i[0], j[0])
|
|---|
| 427 | outfile.write(out_line)
|
|---|
| 428 | break
|
|---|
| 429 |
|
|---|
| 430 | outfile.close()
|
|---|
| 431 |
|
|---|
| 432 | # create reclass of target from reclass rules file
|
|---|
| 433 | result = grass.core.find_file(matched, element='cell')
|
|---|
| 434 | if result['fullname']:
|
|---|
| 435 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 436 | name=matched)
|
|---|
| 437 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 438 | rules=outfile.name)
|
|---|
| 439 | else:
|
|---|
| 440 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 441 | rules=outfile.name)
|
|---|
| 442 |
|
|---|
| 443 | # Cleanup
|
|---|
| 444 | # remove the rules file
|
|---|
| 445 | grass.try_remove(outfile.name)
|
|---|
| 446 |
|
|---|
| 447 | # return reclass of target with histogram that matches original
|
|---|
| 448 | return matched
|
|---|
| 449 |
|
|---|
| 450 | if __name__ == "__main__":
|
|---|
| 451 | options, flags = grass.parser()
|
|---|
| 452 | main()
|
|---|
