← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to scripts/i.pansharpen/i.pansharpen.py

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/i.pansharpen/i.pansharpen.py
 
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", col='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 correspondes 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()