| 1 | #!/usr/bin/env python
|
|---|
| 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 |
|
|---|
| 91 | import os
|
|---|
| 92 |
|
|---|
| 93 | try:
|
|---|
| 94 | import numpy as np
|
|---|
| 95 | hasNumPy = True
|
|---|
| 96 | except ImportError:
|
|---|
| 97 | hasNumPy = False
|
|---|
| 98 |
|
|---|
| 99 | import grass.script as grass
|
|---|
| 100 |
|
|---|
| 101 | # i18N
|
|---|
| 102 | try:
|
|---|
| 103 | import gettext
|
|---|
| 104 | hasGetText = True
|
|---|
| 105 | gettext.install('grassmods', os.path.join(os.getenv("GISBASE"), 'locale'))
|
|---|
| 106 | except ImportError:
|
|---|
| 107 | hasGetText = False
|
|---|
| 108 |
|
|---|
| 109 | def main():
|
|---|
| 110 | if not hasNumPy:
|
|---|
| 111 | grass.fatal(_("Required dependency NumPy not found. Exiting."))
|
|---|
| 112 |
|
|---|
| 113 | sharpen = options['method'] # sharpening algorithm
|
|---|
| 114 | ms1_orig = options['blue'] # blue channel
|
|---|
| 115 | ms2_orig = options['green'] # green channel
|
|---|
| 116 | ms3_orig = options['red'] # red channel
|
|---|
| 117 | pan_orig = options['pan'] # high res pan channel
|
|---|
| 118 | out = options['output'] # prefix for output RGB maps
|
|---|
| 119 | bits = options['bitdepth'] # bit depth of image channels
|
|---|
| 120 | bladjust = flags['l'] # adjust blue channel
|
|---|
| 121 | sproc = flags['s'] # serial processing
|
|---|
| 122 |
|
|---|
| 123 | # Internationalization
|
|---|
| 124 | if not hasGetText:
|
|---|
| 125 | grass.warning(_("No gettext international language support"))
|
|---|
| 126 |
|
|---|
| 127 | # Checking bit depth
|
|---|
| 128 | bits = float(bits)
|
|---|
| 129 | if bits < 2 or bits > 30:
|
|---|
| 130 | grass.warning(_("Bit depth is outside acceptable range"))
|
|---|
| 131 | return
|
|---|
| 132 |
|
|---|
| 133 | outb = grass.core.find_file('%s_blue' % out)
|
|---|
| 134 | outg = grass.core.find_file('%s_green' % out)
|
|---|
| 135 | outr = grass.core.find_file('%s_red' % out)
|
|---|
| 136 |
|
|---|
| 137 | if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and not grass.overwrite():
|
|---|
| 138 | grass.warning(_('Maps with selected output prefix names already exist.'
|
|---|
| 139 | ' Delete them or use overwrite flag'))
|
|---|
| 140 | return
|
|---|
| 141 |
|
|---|
| 142 | pid = str(os.getpid())
|
|---|
| 143 |
|
|---|
| 144 | # convert input image channels to 8 bit for processing
|
|---|
| 145 | ms1 = 'tmp%s_ms1' % pid
|
|---|
| 146 | ms2 = 'tmp%s_ms2' % pid
|
|---|
| 147 | ms3 = 'tmp%s_ms3' % pid
|
|---|
| 148 | pan = 'tmp%s_pan' % pid
|
|---|
| 149 |
|
|---|
| 150 | if bits == 8:
|
|---|
| 151 | grass.message(_("Using 8bit image channels"))
|
|---|
| 152 | if sproc:
|
|---|
| 153 | # serial processing
|
|---|
| 154 | grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
|
|---|
| 155 | quiet=True, overwrite=True)
|
|---|
| 156 | grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
|
|---|
| 157 | quiet=True, overwrite=True)
|
|---|
| 158 | grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
|
|---|
| 159 | quiet=True, overwrite=True)
|
|---|
| 160 | grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
|
|---|
| 161 | quiet=True, overwrite=True)
|
|---|
| 162 | else:
|
|---|
| 163 | # parallel processing
|
|---|
| 164 | pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
|
|---|
| 165 | quiet=True, overwrite=True)
|
|---|
| 166 | pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
|
|---|
| 167 | quiet=True, overwrite=True)
|
|---|
| 168 | pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
|
|---|
| 169 | quiet=True, overwrite=True)
|
|---|
| 170 | pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
|
|---|
| 171 | quiet=True, overwrite=True)
|
|---|
| 172 |
|
|---|
| 173 | pb.wait()
|
|---|
| 174 | pg.wait()
|
|---|
| 175 | pr.wait()
|
|---|
| 176 | pp.wait()
|
|---|
| 177 |
|
|---|
| 178 | else:
|
|---|
| 179 | grass.message(_("Converting image chanels to 8bit for processing"))
|
|---|
| 180 | maxval = pow(2, bits) - 1
|
|---|
| 181 | if sproc:
|
|---|
| 182 | # serial processing
|
|---|
| 183 | grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
|
|---|
| 184 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 185 | grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
|
|---|
| 186 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 187 | grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
|
|---|
| 188 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 189 | grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
|
|---|
| 190 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 191 |
|
|---|
| 192 | else:
|
|---|
| 193 | # parallel processing
|
|---|
| 194 | pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
|
|---|
| 195 | output=ms1, to='0,255', quiet=True, overwrite=True)
|
|---|
| 196 | pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
|
|---|
| 197 | output=ms2, to='0,255', quiet=True, overwrite=True)
|
|---|
| 198 | pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
|
|---|
| 199 | output=ms3, to='0,255', quiet=True, overwrite=True)
|
|---|
| 200 | pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
|
|---|
| 201 | output=pan, to='0,255', quiet=True, overwrite=True)
|
|---|
| 202 |
|
|---|
| 203 | pb.wait()
|
|---|
| 204 | pg.wait()
|
|---|
| 205 | pr.wait()
|
|---|
| 206 | pp.wait()
|
|---|
| 207 |
|
|---|
| 208 | # get PAN resolution:
|
|---|
| 209 | kv = grass.raster_info(map=pan)
|
|---|
| 210 | nsres = kv['nsres']
|
|---|
| 211 | ewres = kv['ewres']
|
|---|
| 212 | panres = (nsres + ewres) / 2
|
|---|
| 213 |
|
|---|
| 214 | # clone current region
|
|---|
| 215 | grass.use_temp_region()
|
|---|
| 216 | grass.run_command('g.region', res=panres, align=pan)
|
|---|
| 217 |
|
|---|
| 218 | # Select sharpening method
|
|---|
| 219 | grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
|
|---|
| 220 | if sharpen == "brovey":
|
|---|
| 221 | brovey(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 222 | elif sharpen == "ihs":
|
|---|
| 223 | ihs(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 224 | elif sharpen == "pca":
|
|---|
| 225 | pca(pan, ms1, ms2, ms3, out, pid, sproc)
|
|---|
| 226 | # Could add other sharpening algorithms here, e.g. wavelet transformation
|
|---|
| 227 |
|
|---|
| 228 | grass.message(_("Assigning grey equalized color tables to output images..."))
|
|---|
| 229 |
|
|---|
| 230 | # equalized grey scales give best contrast
|
|---|
| 231 | grass.message(_("setting pan-sharpened channels to equalized grey scale"))
|
|---|
| 232 | for ch in ['red', 'green', 'blue']:
|
|---|
| 233 | grass.run_command('r.colors', quiet=True, map="%s_%s" % (out, ch),
|
|---|
| 234 | flags="e", color='grey')
|
|---|
| 235 |
|
|---|
| 236 | # Landsat too blue-ish because panchromatic band less sensitive to blue
|
|---|
| 237 | # light, so output blue channed can be modified
|
|---|
| 238 | if bladjust:
|
|---|
| 239 | grass.message(_("Adjusting blue channel color table..."))
|
|---|
| 240 | rules = grass.tempfile()
|
|---|
| 241 | colors = open(rules, 'w')
|
|---|
| 242 | colors.write('5 0 0 0\n20 200 200 200\n40 230 230 230\n67 255 255 255 \n')
|
|---|
| 243 | colors.close()
|
|---|
| 244 |
|
|---|
| 245 | grass.run_command('r.colors', map="%s_blue" % out, rules=rules, quiet=True)
|
|---|
| 246 | os.remove(rules)
|
|---|
| 247 |
|
|---|
| 248 | # output notice
|
|---|
| 249 | grass.verbose(_("The following pan-sharpened output maps have been generated:"))
|
|---|
| 250 | for ch in ['red', 'green', 'blue']:
|
|---|
| 251 | grass.verbose(_("%s_%s") % (out, ch))
|
|---|
| 252 |
|
|---|
| 253 | grass.verbose(_("To visualize output, run: g.region -p raster=%s_red" % out))
|
|---|
| 254 | grass.verbose(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out)))
|
|---|
| 255 | grass.verbose(_("If desired, combine channels into a single RGB map with 'r.composite'."))
|
|---|
| 256 | grass.verbose(_("Channel colors can be rebalanced using i.colors.enhance."))
|
|---|
| 257 |
|
|---|
| 258 | # write cmd history:
|
|---|
| 259 | for ch in ['red', 'green', 'blue']:
|
|---|
| 260 | grass.raster_history("%s_%s" % (out, ch))
|
|---|
| 261 |
|
|---|
| 262 | # create a group with the three outputs
|
|---|
| 263 | #grass.run_command('i.group', group=out,
|
|---|
| 264 | # input="{n}_red,{n}_blue,{n}_green".format(n=out))
|
|---|
| 265 |
|
|---|
| 266 | # Cleanup
|
|---|
| 267 | grass.message(_("cleaning up temp files"))
|
|---|
| 268 | grass.run_command('g.remove', flags="f", type="raster",
|
|---|
| 269 | pattern="tmp%s*" % pid, quiet=True)
|
|---|
| 270 |
|
|---|
| 271 | def brovey(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 272 | grass.verbose(_("Using Brovey algorithm"))
|
|---|
| 273 |
|
|---|
| 274 | # pan/intensity histogram matching using linear regression
|
|---|
| 275 | grass.message(_("Pan channel/intensity histogram matching using linear regression"))
|
|---|
| 276 | outname = 'tmp%s_pan1' % pid
|
|---|
| 277 | panmatch1 = matchhist(pan, ms1, outname)
|
|---|
| 278 |
|
|---|
| 279 | outname = 'tmp%s_pan2' % pid
|
|---|
| 280 | panmatch2 = matchhist(pan, ms2, outname)
|
|---|
| 281 |
|
|---|
| 282 | outname = 'tmp%s_pan3' % pid
|
|---|
| 283 | panmatch3 = matchhist(pan, ms3, outname)
|
|---|
| 284 |
|
|---|
| 285 | outr = '%s_red' % out
|
|---|
| 286 | outg = '%s_green' % out
|
|---|
| 287 | outb = '%s_blue' % out
|
|---|
| 288 |
|
|---|
| 289 | # calculate brovey transformation
|
|---|
| 290 | grass.message(_("Calculating Brovey transformation..."))
|
|---|
| 291 |
|
|---|
| 292 | if sproc:
|
|---|
| 293 | # serial processing
|
|---|
| 294 | e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
|
|---|
| 295 | "$outr" = 1.0 * "$ms3" * "$panmatch3" / k
|
|---|
| 296 | "$outg" = 1.0 * "$ms2" * "$panmatch2" / k
|
|---|
| 297 | "$outb" = 1.0 * "$ms1" * "$panmatch1" / k'''
|
|---|
| 298 | grass.mapcalc(e, outr=outr, outg=outg, outb=outb,
|
|---|
| 299 | panmatch1=panmatch1, panmatch2=panmatch2,
|
|---|
| 300 | panmatch3=panmatch3, ms1=ms1, ms2=ms2, ms3=ms3,
|
|---|
| 301 | overwrite=True)
|
|---|
| 302 | else:
|
|---|
| 303 | # parallel processing
|
|---|
| 304 | pb = grass.mapcalc_start('%s_blue = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 305 | (out, ms1, panmatch1, ms1, ms2, ms3),
|
|---|
| 306 | overwrite=True)
|
|---|
| 307 | pg = grass.mapcalc_start('%s_green = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 308 | (out, ms2, panmatch2, ms1, ms2, ms3),
|
|---|
| 309 | overwrite=True)
|
|---|
| 310 | pr = grass.mapcalc_start('%s_red = (1.0 * %s * %s) / (%s + %s + %s)' %
|
|---|
| 311 | (out, ms3, panmatch3, ms1, ms2, ms3),
|
|---|
| 312 | overwrite=True)
|
|---|
| 313 |
|
|---|
| 314 | pb.wait()
|
|---|
| 315 | pg.wait()
|
|---|
| 316 | pr.wait()
|
|---|
| 317 |
|
|---|
| 318 | # Cleanup
|
|---|
| 319 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 320 | name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|---|
| 321 |
|
|---|
| 322 | def ihs(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 323 | grass.verbose(_("Using IHS<->RGB algorithm"))
|
|---|
| 324 | # transform RGB channels into IHS color space
|
|---|
| 325 | grass.message(_("Transforming to IHS color space..."))
|
|---|
| 326 | grass.run_command('i.rgb.his', overwrite=True,
|
|---|
| 327 | red=ms3,
|
|---|
| 328 | green=ms2,
|
|---|
| 329 | blue=ms1,
|
|---|
| 330 | hue="tmp%s_hue" % pid,
|
|---|
| 331 | intensity="tmp%s_int" % pid,
|
|---|
| 332 | saturation="tmp%s_sat" % pid)
|
|---|
| 333 |
|
|---|
| 334 | # pan/intensity histogram matching using linear regression
|
|---|
| 335 | target = "tmp%s_int" % pid
|
|---|
| 336 | outname = "tmp%s_pan_int" % pid
|
|---|
| 337 | panmatch = matchhist(pan, target, outname)
|
|---|
| 338 |
|
|---|
| 339 | # substitute pan for intensity channel and transform back to RGB color space
|
|---|
| 340 | grass.message(_("Transforming back to RGB color space and sharpening..."))
|
|---|
| 341 | grass.run_command('i.his.rgb', overwrite=True,
|
|---|
| 342 | hue="tmp%s_hue" % pid,
|
|---|
| 343 | intensity="%s" % panmatch,
|
|---|
| 344 | saturation="tmp%s_sat" % pid,
|
|---|
| 345 | red="%s_red" % out,
|
|---|
| 346 | green="%s_green" % out,
|
|---|
| 347 | blue="%s_blue" % out)
|
|---|
| 348 |
|
|---|
| 349 | # Cleanup
|
|---|
| 350 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 351 | name=panmatch)
|
|---|
| 352 |
|
|---|
| 353 | def pca(pan, ms1, ms2, ms3, out, pid, sproc):
|
|---|
| 354 |
|
|---|
| 355 | grass.verbose(_("Using PCA/inverse PCA algorithm"))
|
|---|
| 356 | grass.message(_("Creating PCA images and calculating eigenvectors..."))
|
|---|
| 357 |
|
|---|
| 358 | # initial PCA with RGB channels
|
|---|
| 359 | pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0',
|
|---|
| 360 | input='%s,%s,%s' % (ms1, ms2, ms3),
|
|---|
| 361 | output='tmp%s.pca' % pid)
|
|---|
| 362 | if len(pca_out) < 1:
|
|---|
| 363 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 364 |
|
|---|
| 365 | b1evect = []
|
|---|
| 366 | b2evect = []
|
|---|
| 367 | b3evect = []
|
|---|
| 368 | for l in pca_out.replace('(', ',').replace(')', ',').splitlines():
|
|---|
| 369 | b1evect.append(float(l.split(',')[1]))
|
|---|
| 370 | b2evect.append(float(l.split(',')[2]))
|
|---|
| 371 | b3evect.append(float(l.split(',')[3]))
|
|---|
| 372 |
|
|---|
| 373 | # inverse PCA with hi res pan channel substituted for principal component 1
|
|---|
| 374 | pca1 = 'tmp%s.pca.1' % pid
|
|---|
| 375 | pca2 = 'tmp%s.pca.2' % pid
|
|---|
| 376 | pca3 = 'tmp%s.pca.3' % pid
|
|---|
| 377 | b1evect1 = b1evect[0]
|
|---|
| 378 | b1evect2 = b1evect[1]
|
|---|
| 379 | b1evect3 = b1evect[2]
|
|---|
| 380 | b2evect1 = b2evect[0]
|
|---|
| 381 | b2evect2 = b2evect[1]
|
|---|
| 382 | b2evect3 = b2evect[2]
|
|---|
| 383 | b3evect1 = b3evect[0]
|
|---|
| 384 | b3evect2 = b3evect[1]
|
|---|
| 385 | b3evect3 = b3evect[2]
|
|---|
| 386 |
|
|---|
| 387 | # Histogram matching
|
|---|
| 388 | outname = 'tmp%s_pan1' % pid
|
|---|
| 389 | panmatch1 = matchhist(pan, ms1, outname)
|
|---|
| 390 |
|
|---|
| 391 | outname = 'tmp%s_pan2' % pid
|
|---|
| 392 | panmatch2 = matchhist(pan, ms2, outname)
|
|---|
| 393 |
|
|---|
| 394 | outname = 'tmp%s_pan3' % pid
|
|---|
| 395 | panmatch3 = matchhist(pan, ms3, outname)
|
|---|
| 396 |
|
|---|
| 397 | grass.message(_("Performing inverse PCA ..."))
|
|---|
| 398 |
|
|---|
| 399 | # Get mean value of each channel
|
|---|
| 400 | stats1 = grass.parse_command("r.univar", map=ms1, flags='g',
|
|---|
| 401 | parse=(grass.parse_key_val,
|
|---|
| 402 | {'sep': '='}))
|
|---|
| 403 | stats2 = grass.parse_command("r.univar", map=ms2, flags='g',
|
|---|
| 404 | parse=(grass.parse_key_val,
|
|---|
| 405 | {'sep': '='}))
|
|---|
| 406 | stats3 = grass.parse_command("r.univar", map=ms3, flags='g',
|
|---|
| 407 | parse=(grass.parse_key_val,
|
|---|
| 408 | {'sep': '='}))
|
|---|
| 409 |
|
|---|
| 410 | b1mean = float(stats1['mean'])
|
|---|
| 411 | b2mean = float(stats2['mean'])
|
|---|
| 412 | b3mean = float(stats3['mean'])
|
|---|
| 413 |
|
|---|
| 414 | if sproc:
|
|---|
| 415 | # serial processing
|
|---|
| 416 | outr = '%s_red' % out
|
|---|
| 417 | outg = '%s_green' % out
|
|---|
| 418 | outb = '%s_blue' % out
|
|---|
| 419 |
|
|---|
| 420 | cmd1 = "$outb = (1.0 * $panmatch1 * $b1evect1) + ($pca2 * $b1evect2) + ($pca3 * $b1evect3) + $b1mean"
|
|---|
| 421 | cmd2 = "$outg = (1.0 * $panmatch2 * $b2evect1) + ($pca2 * $b2evect2) + ($pca3 * $b2evect3) + $b2mean"
|
|---|
| 422 | cmd3 = "$outr = (1.0 * $panmatch3 * $b3evect1) + ($pca2 * $b3evect2) + ($pca3 * $b3evect3) + $b3mean"
|
|---|
| 423 |
|
|---|
| 424 | cmd = '\n'.join([cmd1, cmd2, cmd3])
|
|---|
| 425 |
|
|---|
| 426 | grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr,
|
|---|
| 427 | panmatch1=panmatch1,
|
|---|
| 428 | panmatch2=panmatch2,
|
|---|
| 429 | panmatch3=panmatch3,
|
|---|
| 430 | pca2=pca2,
|
|---|
| 431 | pca3=pca3,
|
|---|
| 432 | b1evect1=b1evect1,
|
|---|
| 433 | b2evect1=b2evect1,
|
|---|
| 434 | b3evect1=b3evect1,
|
|---|
| 435 | b1evect2=b1evect2,
|
|---|
| 436 | b2evect2=b2evect2,
|
|---|
| 437 | b3evect2=b3evect2,
|
|---|
| 438 | b1evect3=b1evect3,
|
|---|
| 439 | b2evect3=b2evect3,
|
|---|
| 440 | b3evect3=b3evect3,
|
|---|
| 441 | b1mean=b1mean,
|
|---|
| 442 | b2mean=b2mean,
|
|---|
| 443 | b3mean=b3mean,
|
|---|
| 444 | overwrite=True)
|
|---|
| 445 | else:
|
|---|
| 446 | # parallel processing
|
|---|
| 447 | pb = grass.mapcalc_start('%s_blue = (%s * %f) + (%s * %f) + (%s * %f) + %f'
|
|---|
| 448 | % (out,
|
|---|
| 449 | panmatch1,
|
|---|
| 450 | b1evect1,
|
|---|
| 451 | pca2,
|
|---|
| 452 | b1evect2,
|
|---|
| 453 | pca3,
|
|---|
| 454 | b1evect3,
|
|---|
| 455 | b1mean),
|
|---|
| 456 | overwrite=True)
|
|---|
| 457 |
|
|---|
| 458 | pg = grass.mapcalc_start('%s_green = (%s * %f) + (%s * %f) + (%s * %f) + %f'
|
|---|
| 459 | % (out,
|
|---|
| 460 | panmatch2,
|
|---|
| 461 | b2evect1,
|
|---|
| 462 | pca2,
|
|---|
| 463 | b2evect2,
|
|---|
| 464 | pca3,
|
|---|
| 465 | b2evect3,
|
|---|
| 466 | b2mean),
|
|---|
| 467 | overwrite=True)
|
|---|
| 468 |
|
|---|
| 469 | pr = grass.mapcalc_start('%s_red = (%s * %f) + (%s * %f) + (%s * ''%f) + %f'
|
|---|
| 470 | % (out,
|
|---|
| 471 | panmatch3,
|
|---|
| 472 | b3evect1,
|
|---|
| 473 | pca2,
|
|---|
| 474 | b3evect2,
|
|---|
| 475 | pca3,
|
|---|
| 476 | b3evect3,
|
|---|
| 477 | b3mean),
|
|---|
| 478 | overwrite=True)
|
|---|
| 479 |
|
|---|
| 480 | pb.wait()
|
|---|
| 481 | pg.wait()
|
|---|
| 482 | pr.wait()
|
|---|
| 483 |
|
|---|
| 484 | # Cleanup
|
|---|
| 485 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 486 | name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|---|
| 487 | grass.run_command('g.remove', flags='f', quiet=True, type="raster",
|
|---|
| 488 | pattern='tmp%s*' % pid)
|
|---|
| 489 |
|
|---|
| 490 |
|
|---|
| 491 | def matchhist(original, target, matched):
|
|---|
| 492 | # pan/intensity histogram matching using numpy arrays
|
|---|
| 493 | grass.message(_("Histogram matching..."))
|
|---|
| 494 |
|
|---|
| 495 | # input images
|
|---|
| 496 | original = original.split('@')[0]
|
|---|
| 497 | target = target.split('@')[0]
|
|---|
| 498 | images = [original, target]
|
|---|
| 499 |
|
|---|
| 500 | # create a dictionary to hold arrays for each image
|
|---|
| 501 | arrays = {}
|
|---|
| 502 |
|
|---|
| 503 | for i in images:
|
|---|
| 504 | # calculate number of cells for each grey value for for each image
|
|---|
| 505 | stats_out = grass.pipe_command('r.stats', flags='cin', input=i,
|
|---|
| 506 | sep=':')
|
|---|
| 507 | stats = stats_out.communicate()[0].split('\n')[:-1]
|
|---|
| 508 | stats_dict = dict(s.split(':', 1) for s in stats)
|
|---|
| 509 | total_cells = 0 # total non-null cells
|
|---|
| 510 | for j in stats_dict:
|
|---|
| 511 | stats_dict[j] = int(stats_dict[j])
|
|---|
| 512 | if j != '*':
|
|---|
| 513 | total_cells += stats_dict[j]
|
|---|
| 514 |
|
|---|
| 515 | if total_cells < 1:
|
|---|
| 516 | grass.fatal(_("Input has no data. Check region settings."))
|
|---|
| 517 |
|
|---|
| 518 | # Make a 2x256 structured array for each image with a
|
|---|
| 519 | # cumulative distribution function (CDF) for each grey value.
|
|---|
| 520 | # Grey value is the integer (i4) and cdf is float (f4).
|
|---|
| 521 |
|
|---|
| 522 | arrays[i] = np.zeros((256, ), dtype=('i4,f4'))
|
|---|
| 523 | cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values
|
|---|
| 524 |
|
|---|
| 525 | for n in range(0, 256):
|
|---|
| 526 | if str(n) in stats_dict:
|
|---|
| 527 | num_cells = stats_dict[str(n)]
|
|---|
| 528 | else:
|
|---|
| 529 | num_cells = 0
|
|---|
| 530 |
|
|---|
| 531 | cum_cells += num_cells
|
|---|
| 532 |
|
|---|
| 533 | # cdf is the the number of cells at or below a given grey value
|
|---|
| 534 | # divided by the total number of cells
|
|---|
| 535 | cdf = float(cum_cells) / float(total_cells)
|
|---|
| 536 |
|
|---|
| 537 | # insert values into array
|
|---|
| 538 | arrays[i][n] = (n, cdf)
|
|---|
| 539 |
|
|---|
| 540 | # open file for reclass rules
|
|---|
| 541 | outfile = open(grass.tempfile(), 'w')
|
|---|
| 542 |
|
|---|
| 543 | for i in arrays[original]:
|
|---|
| 544 | # for each grey value and corresponding cdf value in original, find the
|
|---|
| 545 | # cdf value in target that is closest to the target cdf value
|
|---|
| 546 | difference_list = []
|
|---|
| 547 | for j in arrays[target]:
|
|---|
| 548 | # make a list of the difference between each original cdf value and
|
|---|
| 549 | # the target cdf value
|
|---|
| 550 | difference_list.append(abs(i[1] - j[1]))
|
|---|
| 551 |
|
|---|
| 552 | # get the smallest difference in the list
|
|---|
| 553 | min_difference = min(difference_list)
|
|---|
| 554 |
|
|---|
| 555 | for j in arrays[target]:
|
|---|
| 556 | # find the grey value in target that corresponds to the cdf
|
|---|
| 557 | # closest to the original cdf
|
|---|
| 558 | if j[1] == i[1] + min_difference or j[1] == i[1] - min_difference:
|
|---|
| 559 | # build a reclass rules file from the original grey value and
|
|---|
| 560 | # corresponding grey value from target
|
|---|
| 561 | out_line = "%d = %d\n" % (i[0], j[0])
|
|---|
| 562 | outfile.write(out_line)
|
|---|
| 563 | break
|
|---|
| 564 |
|
|---|
| 565 | outfile.close()
|
|---|
| 566 |
|
|---|
| 567 | # create reclass of target from reclass rules file
|
|---|
| 568 | result = grass.core.find_file(matched, element='cell')
|
|---|
| 569 | if result['fullname']:
|
|---|
| 570 | grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|---|
| 571 | name=matched)
|
|---|
| 572 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 573 | rules=outfile.name)
|
|---|
| 574 | else:
|
|---|
| 575 | grass.run_command('r.reclass', input=original, out=matched,
|
|---|
| 576 | rules=outfile.name)
|
|---|
| 577 |
|
|---|
| 578 | # Cleanup
|
|---|
| 579 | # remove the rules file
|
|---|
| 580 | grass.try_remove(outfile.name)
|
|---|
| 581 |
|
|---|
| 582 | # return reclass of target with histogram that matches original
|
|---|
| 583 | return matched
|
|---|
| 584 |
|
|---|
| 585 | if __name__ == "__main__":
|
|---|
| 586 | options, flags = grass.parser()
|
|---|
| 587 | main()
|
|---|
