[GRASS-SVN] r74385 - grass/branches/releasebranch_7_6/scripts/i.pansharpen
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Apr 16 06:38:10 PDT 2019
Author: cmbarton
Date: 2019-04-16 06:38:09 -0700 (Tue, 16 Apr 2019)
New Revision: 74385
Added:
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_brovey542.jpg
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_ihs542.jpg
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_landsat542.jpg
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_pca542.jpg
Modified:
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.html
grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.py
Log:
updating to new version of i.pansharpen merged from trunk
Modified: grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.html
===================================================================
--- grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.html 2019-04-16 13:33:09 UTC (rev 74384)
+++ grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.html 2019-04-16 13:38:09 UTC (rev 74385)
@@ -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,62 +44,71 @@
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>
<h3>Pan sharpening comparison example</h3>
-Pan sharpening of a Landsat image from Boulder, Colorado, USA:
+Pan sharpening of a Landsat image from Boulder, Colorado, USA
+(LANDSAT ETM+ [Landsat 7] spectral bands 5,4,2, and pan band 8):
<div class="code"><pre>
-# 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
- r=p034r032_7dt20010924_z13_30
+# R, G, B composite at 30m
+g.region raster=p034r032_7dt20010924_z13_20 -p
+d.rgb b=p034r032_7dt20010924_z13_20 g=lp034r032_7dt20010924_z13_40
+ r=p034r032_7dt20010924_z13_50
# 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_50 green=p034r032_7dt20010924_z13_40
+ blue=p034r032_7dt20010924_z13_20 pan=p034r032_7dp20010924_z13_80
output=ihs321 method=ihs
# ... likewise with method=brovey and method=pca
# display at 15m
-g.region raster=ihs321_blue -p
-d.rgb b=ihs321_blue g=ihs321_green r=ihs321_red
+g.region raster=ihs542_blue -p
+d.rgb b=ihs542_blue g=ihs542_green r=ihs542_red
</pre></div>
<p>
@@ -109,7 +118,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 +125,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 +134,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 +141,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 +184,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 +201,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 +230,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/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.py
===================================================================
--- grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.py 2019-04-16 13:33:09 UTC (rev 74384)
+++ grass/branches/releasebranch_7_6/scripts/i.pansharpen/i.pansharpen.py 2019-04-16 13:38:09 UTC (rev 74385)
@@ -79,6 +79,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
@@ -107,6 +111,8 @@
out = options['output'] # prefix for output RGB maps
bladjust = flags['l'] # adjust blue channel
sproc = flags['s'] # serial processing
+ rescale = flags['r'] # rescale to spread pixel values to entire 0-255 range
+ rescale = flags['r'] # rescale to spread pixel values to entire 0-255 range
outb = grass.core.find_file('%s_blue' % out)
outg = grass.core.find_file('%s_green' % out)
@@ -125,152 +131,88 @@
ewres = kv['ewres']
panres = (nsres + ewres) / 2
- # clone current region
- grass.use_temp_region()
+ 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)
+ 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),
+ 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)
- grass.run_command('g.region', res=panres, align=pan)
+ pb.wait()
+ pg.wait()
+ pr.wait()
+ pp.wait()
- grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
-
- if sharpen == "brovey":
- grass.verbose(_("Using Brovey algorithm"))
-
- # pan/intensity histogram matching using linear regression
- outname = 'tmp%s_pan1' % pid
- panmatch1 = matchhist(pan, ms1, outname)
-
- outname = 'tmp%s_pan2' % pid
- panmatch2 = matchhist(pan, ms2, outname)
-
- outname = 'tmp%s_pan3' % pid
- panmatch3 = matchhist(pan, ms3, outname)
-
- outr = '%s_red' % out
- outg = '%s_green' % out
- outb = '%s_blue' % out
-
- # calculate brovey transformation
- grass.message(_("Calculating Brovey transformation..."))
-
- if sproc:
- # serial processing
- e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
- "$outr" = 1.0 * "$ms3" * "$panmatch3" / k
- "$outg" = 1.0 * "$ms2" * "$panmatch2" / k
- "$outb" = 1.0 * "$ms1" * "$panmatch1" / k'''
- grass.mapcalc(e, outr=outr, outg=outg, outb=outb,
- panmatch1=panmatch1, panmatch2=panmatch2,
- panmatch3=panmatch3, ms1=ms1, ms2=ms2, ms3=ms3,
- overwrite=True)
else:
- # parallel processing
- pb = grass.mapcalc_start('%s_blue = (1.0 * %s * %s) / (%s + %s + %s)' %
- (out, ms1, panmatch1, ms1, ms2, ms3),
- overwrite=True)
- pg = grass.mapcalc_start('%s_green = (1.0 * %s * %s) / (%s + %s + %s)' %
- (out, ms2, panmatch2, ms1, ms2, ms3),
- overwrite=True)
- pr = grass.mapcalc_start('%s_red = (1.0 * %s * %s) / (%s + %s + %s)' %
- (out, ms3, panmatch3, ms1, ms2, ms3),
- overwrite=True)
+ 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)
- pb.wait()
- pg.wait()
- pr.wait()
+ 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)
- # Cleanup
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
- name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
+ pb.wait()
+ pg.wait()
+ pr.wait()
+ pp.wait()
- elif sharpen == "ihs":
- grass.verbose(_("Using IHS<->RGB algorithm"))
- # transform RGB channels into IHS color space
- grass.message(_("Transforming to IHS color space..."))
- grass.run_command('i.rgb.his', overwrite=True,
- red=ms3,
- green=ms2,
- blue=ms1,
- hue="tmp%s_hue" % pid,
- intensity="tmp%s_int" % pid,
- saturation="tmp%s_sat" % pid)
+ else:
+ grass.message(_("Rescaling image chanels to 8bit for processing"))
- # pan/intensity histogram matching using linear regression
- target = "tmp%s_int" % pid
- outname = "tmp%s_pan_int" % pid
- panmatch = matchhist(pan, target, outname)
+ 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'])
- # substitute pan for intensity channel and transform back to RGB color space
- grass.message(_("Transforming back to RGB color space and sharpening..."))
- grass.run_command('i.his.rgb', overwrite=True,
- hue="tmp%s_hue" % pid,
- intensity="%s" % panmatch,
- saturation="tmp%s_sat" % pid,
- red="%s_red" % out,
- green="%s_green" % out,
- blue="%s_blue" % out)
-
- # Cleanup
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
- name=panmatch)
-
- elif sharpen == "pca":
- grass.verbose(_("Using PCA/inverse PCA algorithm"))
- grass.message(_("Creating PCA images and calculating eigenvectors..."))
-
- # initial PCA with RGB channels
- pca_out = grass.read_command('i.pca', quiet=True, rescale='0,0',
- input='%s,%s,%s' % (ms1, ms2, ms3),
- output='tmp%s.pca' % pid)
- if len(pca_out) < 1:
- grass.fatal(_("Input has no data. Check region settings."))
-
- b1evect = []
- b2evect = []
- b3evect = []
- for l in pca_out.replace('(', ',').replace(')', ',').splitlines():
- b1evect.append(float(l.split(',')[1]))
- b2evect.append(float(l.split(',')[2]))
- b3evect.append(float(l.split(',')[3]))
-
- # inverse PCA with hi res pan channel substituted for principal component 1
- pca1 = 'tmp%s.pca.1' % pid
- pca2 = 'tmp%s.pca.2' % pid
- pca3 = 'tmp%s.pca.3' % pid
- b1evect1 = b1evect[0]
- b1evect2 = b1evect[1]
- b1evect3 = b1evect[2]
- b2evect1 = b2evect[0]
- b2evect2 = b2evect[1]
- b2evect3 = b2evect[2]
- b3evect1 = b3evect[0]
- b3evect2 = b3evect[1]
- b3evect3 = b3evect[2]
-
- outname = 'tmp%s_pan' % pid
- panmatch = matchhist(pan, ms1, outname)
-
- grass.message(_("Performing inverse PCA ..."))
-
- stats1 = grass.parse_command("r.univar", map=ms1, flags='g',
- parse=(grass.parse_key_val,
- {'sep': '='}))
- stats2 = grass.parse_command("r.univar", map=ms2, flags='g',
- parse=(grass.parse_key_val,
- {'sep': '='}))
- stats3 = grass.parse_command("r.univar", map=ms3, flags='g',
- parse=(grass.parse_key_val,
- {'sep': '='}))
-
- b1mean = float(stats1['mean'])
- b2mean = float(stats2['mean'])
- b3mean = float(stats3['mean'])
-
+ maxval = pow(2, bits) - 1
if sproc:
# serial processing
- e = '''eval(k = "$ms1" + "$ms2" + "$ms3")
- "$outr" = 1.0 * "$ms3" * "$panmatch3" / k
- "$outg" = 1.0 * "$ms2" * "$panmatch2" / k
- "$outb" = 1.0* "$ms1" * "$panmatch1" / k'''
+ 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_='%f,%f' % (min_ms2, max_ms2),
+ output=ms2, to='0,255', quiet=True, overwrite=True)
+ 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_='%f,%f' % (min_pan, max_pan),
+ output=pan, to='0,255', quiet=True, overwrite=True)
outr = '%s_red' % out
outg = '%s_green' % out
@@ -291,10 +233,14 @@
overwrite=True)
else:
# parallel processing
- pb = grass.mapcalc_start('%s_blue = (%s * %f) + (%s * %f) + (%s * %f) + %f'
- % (out, panmatch, b1evect1, pca2,
- b2evect1, pca3, b3evect1, b1mean),
- overwrite=True)
+ 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_='%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_='%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_='%f,%f' % (min_pan, max_pan),
+ output=pan, to='0,255', quiet=True, overwrite=True)
pg = grass.mapcalc_start('%s_green = (%s * %f) + (%s * %f) + (%s * %f) + %f'
% (out, panmatch, b1evect2, pca2,
@@ -301,11 +247,13 @@
b2evect2, pca3, b3evect2, b2mean),
overwrite=True)
- pr = grass.mapcalc_start('%s_red = (%s * %f) + (%s * %f) + (%s * ''%f) + %f'
- % (out, panmatch, b1evect3, pca2,
- b2evect3, pca3, b3evect3, b3mean),
- overwrite=True)
+ # get PAN resolution:
+ kv = grass.raster_info(map=pan)
+ nsres = kv['nsres']
+ ewres = kv['ewres']
+ panres = (nsres + ewres) / 2
+
pr.wait()
pg.wait()
pb.wait()
@@ -326,14 +274,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)
- os.remove(rules)
-
# output notice
grass.verbose(_("The following pan-sharpened output maps have been generated:"))
for ch in ['red', 'green', 'blue']:
@@ -369,9 +316,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)
@@ -388,7 +335,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):
@@ -404,7 +351,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')
Copied: grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_brovey542.jpg (from rev 74382, grass/trunk/scripts/i.pansharpen/i_pansharpen_rgb_brovey542.jpg)
===================================================================
(Binary files differ)
Copied: grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_ihs542.jpg (from rev 74382, grass/trunk/scripts/i.pansharpen/i_pansharpen_rgb_ihs542.jpg)
===================================================================
(Binary files differ)
Copied: grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_landsat542.jpg (from rev 74382, grass/trunk/scripts/i.pansharpen/i_pansharpen_rgb_landsat542.jpg)
===================================================================
(Binary files differ)
Copied: grass/branches/releasebranch_7_6/scripts/i.pansharpen/i_pansharpen_rgb_pca542.jpg (from rev 74382, grass/trunk/scripts/i.pansharpen/i_pansharpen_rgb_pca542.jpg)
===================================================================
(Binary files differ)
More information about the grass-commit
mailing list