[GRASS-SVN] r74381 - grass/trunk/scripts/i.pansharpen
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 16 02:30:36 PDT 2019
Author: cmbarton
Date: 2019-04-16 02:30:36 -0700 (Tue, 16 Apr 2019)
New Revision: 74381
Modified:
grass/trunk/scripts/i.pansharpen/i.pansharpen.html
grass/trunk/scripts/i.pansharpen/i.pansharpen.py
Log:
added new channel stretch option and updated manual
Modified: grass/trunk/scripts/i.pansharpen/i.pansharpen.html
===================================================================
--- grass/trunk/scripts/i.pansharpen/i.pansharpen.html 2019-04-16 01:12:49 UTC (rev 74380)
+++ grass/trunk/scripts/i.pansharpen/i.pansharpen.html 2019-04-16 09:30:36 UTC (rev 74381)
@@ -5,39 +5,39 @@
lower resolution bands can then be combined into an RGB color image at a
higher (more detailed) resolution than is possible using the original 3
bands. For example, Landsat ETM has low resolution spectral bands 1 (blue),
-2 (green), 3 (red), 4 (near IR), 5 (mid-IR), and 7 (mid-IR) at 30m resolution,
+2 (green), 3 (red), 4 (near IR), 5 (mid-IR), and 7 (mid-IR) at 30m resolution,
and a high resolution panchromatic band 8 at 15m resolution. Pan sharpening
-allows bands 3-2-1 (or other combinations of 30m resolution bands like 4-3-2
+allows bands 3-2-1 (or other combinations of 30m resolution bands like 4-3-2
or 5-4-2) to be combined into a 15m resolution color image.
<br><br>
-i.pansharpen offers a choice of three different 'pan sharpening'
+i.pansharpen offers a choice of three different 'pan sharpening'
algorithms: IHS, Brovey, and PCA.
<br><br>
-For <em>IHS pan sharpening</em>, the original 3 lower resolution bands, selected
-as red, green and blue channels for creating an RGB composite image, are
-transformed into IHS (intensity, hue, and saturation) color space. The
-panchromatic band is then substituted for the intensity channel (I), combined
-with the original hue (H) and saturation (S) channels, and transformed back to
-RGB color space at the higher resolution of the panchromatic band. The
+For <em>IHS pan sharpening</em>, the original 3 lower resolution bands, selected
+as red, green and blue channels for creating an RGB composite image, are
+transformed into IHS (intensity, hue, and saturation) color space. The
+panchromatic band is then substituted for the intensity channel (I), combined
+with the original hue (H) and saturation (S) channels, and transformed back to
+RGB color space at the higher resolution of the panchromatic band. The
algorithm for this can be represented as: RGB -> IHS -> [pan]HS -> RGB.
<br><br>
-With a <em>Brovey pan sharpening</em>, each of the 3 lower resolution bands and
-panchromatic band are combined using the following algorithm to calculate
+With a <em>Brovey pan sharpening</em>, each of the 3 lower resolution bands and
+panchromatic band are combined using the following algorithm to calculate
3 new bands at the higher resolution (example for band 1):
<pre>
- band1
+ band1
new band1 = ----------------------- * panband
band1 + band2 + band3
</pre>
-In <em>PCA pan sharpening</em>, a principal component analysis is performed on the
+In <em>PCA pan sharpening</em>, a principal component analysis is performed on the
original 3 lower resolution bands to create 3 principal component images
(PC1, PC2, and PC3) and their associated eigenvectors (EV), such that:
<pre>
-
+
band1 band2 band3
PC1: EV1-1 EV1-2 EV1-3
PC2: EV2-1 EV2-2 EV2-3
-PC3: EV3-1 EV3-2 EV3-3
+PC3: EV3-1 EV3-2 EV3-3
and
@@ -44,39 +44,47 @@
PC1 = EV1-1 * band1 + EV1-2 * band2 + EV1-3 * band3 - mean(bands 1,2,3)
</pre>
-An inverse PCA is then performed, substituting the panchromatic band for PC1.
-To do this, the eigenvectors matrix is inverted (in this case transposed), the
-PC images are multiplied by the eigenvectors with the panchromatic band
-substituted for PC1, and mean of each band is added to each transformed image
+An inverse PCA is then performed, substituting the panchromatic band for PC1.
+To do this, the eigenvectors matrix is inverted (in this case transposed), the
+PC images are multiplied by the eigenvectors with the panchromatic band
+substituted for PC1, and mean of each band is added to each transformed image
band using the following algorithm (example for band 1):
<pre>
-band1' = pan * EV1-1 + PC2 * EV2-1 + PC3 * EV3-1 + mean(band1)
-
+band1 = pan * EV1-1 + PC2 * EV1-2 + PC3 * EV1-3 + mean(band1)
+
</pre>
-The assignment of the channels depends on the satellite. Examples of satellite
-imagery with high resolution panchromatic bands, and lower resolution spectral
+The assignment of the channels depends on the satellite. Examples of satellite
+imagery with high resolution panchromatic bands, and lower resolution spectral
bands include Landsat 7 ETM, QuickBird, and SPOT.
<br>
<h2>NOTES</h2>
-The module currently only works for 8-bit images.
+The module works for 2-bit to 30-bit images. All images are rescaled to 8-bit
+for processing. By default, the entire possible range for the selected bit depth is
+rescaled to 8-bit. For example, the range of 0-65535 for a 16-bit image is
+rescaled to 0-255). The 'r' flag allows the range of pixel values actually
+present in an image rescaled to a full 8-bit range. For example, a 16 bit image
+might only have pixels that range from 70 to 35000; this range of 70-35000 would
+be rescaled to 0-255. This can give better visual distinction to features,
+especially when the range of actual values in an image only occupies a
+relatively limited portion of the possible range.
<br><br>
-The command temporarily changes the computational region to the high
-resolution of the panchromatic band during sharpening calculations, then
-restores the previous region settings. The current region coordinates (and
-null values) are respected. The high resolution panchromatic image is
-histogram matched to the band it is replaces prior to substitution (i.e., the
-intensity channel for IHS sharpening, the low res band selected for each color
+i.pansharpen temporarily changes the computational region to the high
+resolution of the panchromatic band during sharpening calculations, then
+restores the previous region settings. The current region coordinates (and
+null values) are respected. The high resolution panchromatic image is
+histogram matched to the band it is replaces prior to substitution (i.e., the
+intensity channel for IHS sharpening, the low res band selected for each color
channel with Brovey sharpening, and the PC1 image for PCA sharpening).
<br><br>
-By default, the command will attempt to employ parallel processing, using
-up to 3 cores simultaneously. The -s flag will disable parallel processing,
+By default, the command will attempt to employ parallel processing, using
+up to 3 cores simultaneously. The -s flag will disable parallel processing,
but does use an optimized r.mapcalc expression to reduce disk I/O.
<br><br>
-The three pan-sharpened output channels may be combined with <em>d.rgb</em> or
+The three pan-sharpened output channels may be combined with <em>d.rgb</em> or
<em>r.composite</em>. Colors may be optionally optimized with <em>i.colors.enhance</em>.
-While the resulting color image will be at the higher resolution in all cases,
-the 3 pan sharpening algorithms differ in terms of spectral response.
+While the resulting color image will be at the higher resolution in all cases,
+the 3 pan sharpening algorithms differ in terms of spectral response.
<h2>EXAMPLES</h2>
@@ -85,14 +93,14 @@
Pan sharpening of a Landsat image from Boulder, Colorado, USA:
<div class="code"><pre>
-# R, G, B composite at 30m
+# R, G, B composite at 30m
g.region raster=p034r032_7dt20010924_z13_10 -p
-d.rgb b=p034r032_7dt20010924_z13_10 g=lp034r032_7dt20010924_z13_20
+d.rgb b=p034r032_7dt20010924_z13_10 g=lp034r032_7dt20010924_z13_20
r=p034r032_7dt20010924_z13_30
# i.pansharpen with IHS algorithm
-i.pansharpen red=p034r032_7dt20010924_z13_30 green=p034r032_7dt20010924_z13_20
- blue=p034r032_7dt20010924_z13_10 pan=p034r032_7dp20010924_z13_80
+i.pansharpen red=p034r032_7dt20010924_z13_30 green=p034r032_7dt20010924_z13_20
+ blue=p034r032_7dt20010924_z13_10 pan=p034r032_7dp20010924_z13_80
output=ihs321 method=ihs
# ... likewise with method=brovey and method=pca
@@ -109,7 +117,7 @@
<table border=1>
<tr>
<td align=center>
- <img src="i_pansharpen_rgb_landsat321.jpg" alt="R, G, B composite of Landsat at 30m">
+ <img src="i_pansharpen_rgb_landsat542.jpg" alt="R, G, B composite of Landsat at 30m">
<br>
<font size="-1">
<i>R, G, B composite of Landsat at 30m</i>
@@ -116,7 +124,7 @@
</font>
</td>
<td align=center>
- <img src="i_pansharpen_rgb_brovey321.jpg" alt="R, G, B composite of Brovey sharpened image at 15m">
+ <img src="i_pansharpen_rgb_brovey542.jpg" alt="R, G, B composite of Brovey sharpened image at 15m">
<br>
<font size="-1">
<i>R, G, B composite of Brovey sharpened image at 15m</i>
@@ -125,7 +133,7 @@
</tr>
<tr>
<td align=center>
- <img src="i_pansharpen_rgb_ihs321.jpg" alt="R, G, B composite of IHS sharpened image at 15m">
+ <img src="i_pansharpen_rgb_ihs542.jpg" alt="R, G, B composite of IHS sharpened image at 15m">
<br>
<font size="-1">
<i>R, G, B composite of IHS sharpened image at 15m</i>
@@ -132,7 +140,7 @@
</font>
</td>
<td align=center>
- <img src="i_pansharpen_rgb_pca321.jpg" alt="R, G, B composite of PCA sharpened image at 15m">
+ <img src="i_pansharpen_rgb_pca542.jpg" alt="R, G, B composite of PCA sharpened image at 15m">
<br>
<font size="-1">
<i>R, G, B composite of PCA sharpened image at 15m"</i>
@@ -175,7 +183,6 @@
<h2>SEE ALSO</h2>
<em>
-<a href="i.colors.enhance.html">i.colors.enhance</a>,
<a href="i.his.rgb.html">i.his.rgb</a>,
<a href="i.rgb.his.html">i.rgb.his</a>,
<a href="i.pca.html">i.pca</a>,
@@ -193,26 +200,26 @@
Proc. of the 14th International Symposium on Remote Sensing
of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007
-<li>Amarsaikhan, D., Douglas, T. (2004). Data fusion and multisource image
+<li>Amarsaikhan, D., Douglas, T. (2004). Data fusion and multisource image
classification. International Journal of Remote Sensing, 25(17), 3529-3539.
-<li>Behnia, P. (2005). Comparison between four methods for data fusion of ETM+
+<li>Behnia, P. (2005). Comparison between four methods for data fusion of ETM+
multispectral and pan images. Geo-spatial Information Science, 8(2), 98-103.
-
-<li>Du, Q., Younan, N. H., King, R., Shah, V. P. (2007). On the Performance
- Evaluation of Pan-Sharpening Techniques. Geoscience and Remote Sensing
+
+<li>Du, Q., Younan, N. H., King, R., Shah, V. P. (2007). On the Performance
+ Evaluation of Pan-Sharpening Techniques. Geoscience and Remote Sensing
Letters, IEEE, 4(4), 518-522.
-<li>Karathanassi, V., Kolokousis, P., Ioannidou, S. (2007). A comparison
- study on fusion methods using evaluation indicators. International Journal
+<li>Karathanassi, V., Kolokousis, P., Ioannidou, S. (2007). A comparison
+ study on fusion methods using evaluation indicators. International Journal
of Remote Sensing, 28(10), 2309-2341.
<li>Neteler, M, D. Grasso, I. Michelazzi, L. Miori, S. Merler, and C.
- Furlanello (2005). An integrated toolbox for image registration, fusion and
+ Furlanello (2005). An integrated toolbox for image registration, fusion and
classification. International Journal of Geoinformatics, 1(1):51-61
(<a href="http://www.grassbook.org/wp-content/uploads/neteler/papers/neteler2005_IJG_051-061_draft.pdf">PDF</a>)
-<li>Pohl, C, and J.L van Genderen (1998). Multisensor image fusion in remote
+<li>Pohl, C, and J.L van Genderen (1998). Multisensor image fusion in remote
sensing: concepts, methods and application. Int. J. of Rem. Sens., 19, 823-854.
</ul>
@@ -222,6 +229,6 @@
Michael Barton (Arizona State University, USA)<br>
with contributions from Markus Neteler (ITC-irst, Italy); Glynn Clements;
-Luca Delucchi (Fondazione E. Mach, Italy); Markus Metz; and Hamish Bowman.
+Luca Delucchi (Fondazione E. Mach, Italy); Markus Metz; and Hamish Bowman.
<p><i>Last changed: $Date$</i>
Modified: grass/trunk/scripts/i.pansharpen/i.pansharpen.py
===================================================================
--- grass/trunk/scripts/i.pansharpen/i.pansharpen.py 2019-04-16 01:12:49 UTC (rev 74380)
+++ grass/trunk/scripts/i.pansharpen/i.pansharpen.py 2019-04-16 09:30:36 UTC (rev 74381)
@@ -87,6 +87,10 @@
#% key: l
#% description: Rebalance blue channel for LANDSAT
#%end
+#%flag
+#% key: r
+#% description: Rescale (stretch) the range of pixel values in each channel to the entire 0-255 8-bit range for processing (see notes)
+#%end
import os
@@ -112,6 +116,7 @@
bits = options['bitdepth'] # bit depth of image channels
bladjust = flags['l'] # adjust blue channel
sproc = flags['s'] # serial processing
+ rescale = flags['r'] # rescale to spread pixel values to entire 0-255 range
# Checking bit depth
bits = float(bits)
@@ -136,57 +141,98 @@
ms3 = 'tmp%s_ms3' % pid
pan = 'tmp%s_pan' % pid
- if bits == 8:
- grass.message(_("Using 8bit image channels"))
- if sproc:
- # serial processing
- grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
- quiet=True, overwrite=True)
- grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
- quiet=True, overwrite=True)
- grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
- quiet=True, overwrite=True)
- grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
- quiet=True, overwrite=True)
- else:
- # parallel processing
- pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
+ if rescale == False:
+ if bits == 8:
+ grass.message(_("Using 8bit image channels"))
+ if sproc:
+ # serial processing
+ grass.run_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
quiet=True, overwrite=True)
- pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
+ grass.run_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
quiet=True, overwrite=True)
- pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
+ grass.run_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
quiet=True, overwrite=True)
- pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
+ grass.run_command('g.copy', raster='%s,%s' % (pan_orig, pan),
quiet=True, overwrite=True)
+ else:
+ # parallel processing
+ pb = grass.start_command('g.copy', raster='%s,%s' % (ms1_orig, ms1),
+ quiet=True, overwrite=True)
+ pg = grass.start_command('g.copy', raster='%s,%s' % (ms2_orig, ms2),
+ quiet=True, overwrite=True)
+ pr = grass.start_command('g.copy', raster='%s,%s' % (ms3_orig, ms3),
+ quiet=True, overwrite=True)
+ pp = grass.start_command('g.copy', raster='%s,%s' % (pan_orig, pan),
+ quiet=True, overwrite=True)
- pb.wait()
- pg.wait()
- pr.wait()
- pp.wait()
+ pb.wait()
+ pg.wait()
+ pr.wait()
+ pp.wait()
+ else:
+ grass.message(_("Converting image chanels to 8bit for processing"))
+ maxval = pow(2, bits) - 1
+ if sproc:
+ # serial processing
+ grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+ output=ms1, to='0,255', quiet=True, overwrite=True)
+ grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+ output=ms2, to='0,255', quiet=True, overwrite=True)
+ grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+ output=ms3, to='0,255', quiet=True, overwrite=True)
+ grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+ output=pan, to='0,255', quiet=True, overwrite=True)
+
+ else:
+ # parallel processing
+ pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+ output=ms1, to='0,255', quiet=True, overwrite=True)
+ pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+ output=ms2, to='0,255', quiet=True, overwrite=True)
+ pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+ output=ms3, to='0,255', quiet=True, overwrite=True)
+ pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+ output=pan, to='0,255', quiet=True, overwrite=True)
+
+ pb.wait()
+ pg.wait()
+ pr.wait()
+ pp.wait()
+
else:
- grass.message(_("Converting image chanels to 8bit for processing"))
+ grass.message(_("Rescaling image chanels to 8bit for processing"))
+
+ min_ms1 = int(grass.raster_info(ms1_orig)['min'])
+ max_ms1 = int(grass.raster_info(ms1_orig)['max'])
+ min_ms2 = int(grass.raster_info(ms2_orig)['min'])
+ max_ms2 = int(grass.raster_info(ms2_orig)['max'])
+ min_ms3 = int(grass.raster_info(ms3_orig)['min'])
+ max_ms3 = int(grass.raster_info(ms3_orig)['max'])
+ min_pan = int(grass.raster_info(pan_orig)['min'])
+ max_pan = int(grass.raster_info(pan_orig)['max'])
+
maxval = pow(2, bits) - 1
if sproc:
# serial processing
- grass.run_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+ grass.run_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
output=ms1, to='0,255', quiet=True, overwrite=True)
- grass.run_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+ grass.run_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
output=ms2, to='0,255', quiet=True, overwrite=True)
- grass.run_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+ grass.run_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
output=ms3, to='0,255', quiet=True, overwrite=True)
- grass.run_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+ grass.run_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
output=pan, to='0,255', quiet=True, overwrite=True)
else:
# parallel processing
- pb = grass.start_command('r.rescale', input=ms1_orig, from_='0,%f' % maxval,
+ pb = grass.start_command('r.rescale', input=ms1_orig, from_='%f,%f' % (min_ms1, max_ms1),
output=ms1, to='0,255', quiet=True, overwrite=True)
- pg = grass.start_command('r.rescale', input=ms2_orig, from_='0,%f' % maxval,
+ pg = grass.start_command('r.rescale', input=ms2_orig, from_='%f,%f' % (min_ms2, max_ms2),
output=ms2, to='0,255', quiet=True, overwrite=True)
- pr = grass.start_command('r.rescale', input=ms3_orig, from_='0,%f' % maxval,
+ pr = grass.start_command('r.rescale', input=ms3_orig, from_='%f,%f' % (min_ms3, max_ms3),
output=ms3, to='0,255', quiet=True, overwrite=True)
- pp = grass.start_command('r.rescale', input=pan_orig, from_='0,%f' % maxval,
+ pp = grass.start_command('r.rescale', input=pan_orig, from_='%f,%f' % (min_pan, max_pan),
output=pan, to='0,255', quiet=True, overwrite=True)
pb.wait()
@@ -194,6 +240,7 @@
pr.wait()
pp.wait()
+
# get PAN resolution:
kv = grass.raster_info(map=pan)
nsres = kv['nsres']
@@ -226,14 +273,13 @@
# light, so output blue channed can be modified
if bladjust:
grass.message(_("Adjusting blue channel color table..."))
- rules = grass.tempfile()
- colors = open(rules, 'w')
- colors.write('5 0 0 0\n20 200 200 200\n40 230 230 230\n67 255 255 255 \n')
- colors.close()
+ blue_colors = ['0 0 0 0\n5% 0 0 0\n67% 255 255 255\n100% 255 255 255']
+ # these previous colors are way too blue for landsat
+ # 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']
+ bc = grass.feed_command('r.colors', quiet = True, map = "%s_blue" % out, rules = "-")
+ bc.stdin.write('\n'.join(blue_colors))
+ bc.stdin.close()
- grass.run_command('r.colors', map="%s_blue" % out, rules=rules, quiet=True)
- os.remove(rules)
-
# output notice
grass.verbose(_("The following pan-sharpened output maps have been generated:"))
for ch in ['red', 'green', 'blue']:
@@ -499,9 +545,9 @@
# create a dictionary to hold arrays for each image
arrays = {}
- for i in images:
+ for img in images:
# calculate number of cells for each grey value for for each image
- stats_out = grass.pipe_command('r.stats', flags='cin', input=i,
+ stats_out = grass.pipe_command('r.stats', flags='cin', input=img,
sep=':')
stats = stats_out.communicate()[0].split('\n')[:-1]
stats_dict = dict(s.split(':', 1) for s in stats)
@@ -518,7 +564,7 @@
# cumulative distribution function (CDF) for each grey value.
# Grey value is the integer (i4) and cdf is float (f4).
- arrays[i] = np.zeros((256, ), dtype=('i4,f4'))
+ arrays[img] = np.zeros((256, ), dtype=('i4,f4'))
cum_cells = 0 # cumulative total of cells for sum of current and all lower grey values
for n in range(0, 256):
@@ -534,7 +580,7 @@
cdf = float(cum_cells) / float(total_cells)
# insert values into array
- arrays[i][n] = (n, cdf)
+ arrays[img][n] = (n, cdf)
# open file for reclass rules
outfile = open(grass.tempfile(), 'w')
More information about the grass-commit
mailing list