|
@@ -16,9 +16,9 @@
|
|
|
#
|
|
|
# COPYRIGHT: (C) 2002-2019 by the GRASS Development Team
|
|
|
#
|
|
|
-# This program is free software under the GNU General Public
|
|
|
-# License (>=v2). Read the file COPYING that comes with GRASS
|
|
|
-# for details.
|
|
|
+# This program is free software under the GNU General Public
|
|
|
+# License (>=v2). Read the file COPYING that comes with GRASS
|
|
|
+# for details.
|
|
|
#
|
|
|
# REFERENCES:
|
|
|
# Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS and merged MSS/RBV
|
|
@@ -96,6 +96,7 @@ import os
|
|
|
|
|
|
try:
|
|
|
import numpy as np
|
|
|
+
|
|
|
hasNumPy = True
|
|
|
except ImportError:
|
|
|
hasNumPy = False
|
|
@@ -107,16 +108,16 @@ def main():
|
|
|
if not hasNumPy:
|
|
|
grass.fatal(_("Required dependency NumPy not found. Exiting."))
|
|
|
|
|
|
- sharpen = options['method'] # sharpening algorithm
|
|
|
- ms1_orig = options['blue'] # blue channel
|
|
|
- ms2_orig = options['green'] # green channel
|
|
|
- ms3_orig = options['red'] # red channel
|
|
|
- pan_orig = options['pan'] # high res pan channel
|
|
|
- out = options['output'] # prefix for output RGB maps
|
|
|
- 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
|
|
|
+ sharpen = options["method"] # sharpening algorithm
|
|
|
+ ms1_orig = options["blue"] # blue channel
|
|
|
+ ms2_orig = options["green"] # green channel
|
|
|
+ ms3_orig = options["red"] # red channel
|
|
|
+ pan_orig = options["pan"] # high res pan channel
|
|
|
+ out = options["output"] # prefix for output RGB maps
|
|
|
+ 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)
|
|
@@ -124,46 +125,84 @@ def main():
|
|
|
grass.warning(_("Bit depth is outside acceptable range"))
|
|
|
return
|
|
|
|
|
|
- outb = grass.core.find_file('%s_blue' % out)
|
|
|
- outg = grass.core.find_file('%s_green' % out)
|
|
|
- outr = grass.core.find_file('%s_red' % out)
|
|
|
-
|
|
|
- if (outb['name'] != '' or outg['name'] != '' or outr['name'] != '') and not grass.overwrite():
|
|
|
- grass.warning(_('Maps with selected output prefix names already exist.'
|
|
|
- ' Delete them or use overwrite flag'))
|
|
|
+ outb = grass.core.find_file("%s_blue" % out)
|
|
|
+ outg = grass.core.find_file("%s_green" % out)
|
|
|
+ outr = grass.core.find_file("%s_red" % out)
|
|
|
+
|
|
|
+ if (
|
|
|
+ outb["name"] != "" or outg["name"] != "" or outr["name"] != ""
|
|
|
+ ) and not grass.overwrite():
|
|
|
+ grass.warning(
|
|
|
+ _(
|
|
|
+ "Maps with selected output prefix names already exist."
|
|
|
+ " Delete them or use overwrite flag"
|
|
|
+ )
|
|
|
+ )
|
|
|
return
|
|
|
|
|
|
pid = str(os.getpid())
|
|
|
|
|
|
# convert input image channels to 8 bit for processing
|
|
|
- ms1 = 'tmp%s_ms1' % pid
|
|
|
- ms2 = 'tmp%s_ms2' % pid
|
|
|
- ms3 = 'tmp%s_ms3' % pid
|
|
|
- pan = 'tmp%s_pan' % pid
|
|
|
+ ms1 = "tmp%s_ms1" % pid
|
|
|
+ ms2 = "tmp%s_ms2" % pid
|
|
|
+ ms3 = "tmp%s_ms3" % pid
|
|
|
+ pan = "tmp%s_pan" % pid
|
|
|
|
|
|
if not rescale:
|
|
|
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)
|
|
|
+ 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)
|
|
|
+ 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()
|
|
@@ -175,25 +214,81 @@ def main():
|
|
|
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)
|
|
|
+ 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 = 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()
|
|
@@ -203,53 +298,108 @@ def main():
|
|
|
else:
|
|
|
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'])
|
|
|
+ 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_='%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)
|
|
|
+ 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,
|
|
|
+ )
|
|
|
|
|
|
else:
|
|
|
# parallel processing
|
|
|
- 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)
|
|
|
+ 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,
|
|
|
+ )
|
|
|
|
|
|
pb.wait()
|
|
|
pg.wait()
|
|
|
pr.wait()
|
|
|
pp.wait()
|
|
|
|
|
|
-
|
|
|
# get PAN resolution:
|
|
|
kv = grass.raster_info(map=pan)
|
|
|
- nsres = kv['nsres']
|
|
|
- ewres = kv['ewres']
|
|
|
+ nsres = kv["nsres"]
|
|
|
+ ewres = kv["ewres"]
|
|
|
panres = (nsres + ewres) / 2
|
|
|
|
|
|
# clone current region
|
|
|
grass.use_temp_region()
|
|
|
- grass.run_command('g.region', res=panres, align=pan)
|
|
|
+ grass.run_command("g.region", res=panres, align=pan)
|
|
|
|
|
|
# Select sharpening method
|
|
|
grass.message(_("Performing pan sharpening with hi res pan image: %f" % panres))
|
|
@@ -265,89 +415,109 @@ def main():
|
|
|
|
|
|
# equalized grey scales give best contrast
|
|
|
grass.message(_("setting pan-sharpened channels to equalized grey scale"))
|
|
|
- for ch in ['red', 'green', 'blue']:
|
|
|
- grass.run_command('r.colors', quiet=True, map="%s_%s" % (out, ch),
|
|
|
- flags="e", color='grey')
|
|
|
+ for ch in ["red", "green", "blue"]:
|
|
|
+ grass.run_command(
|
|
|
+ "r.colors", quiet=True, map="%s_%s" % (out, ch), flags="e", color="grey"
|
|
|
+ )
|
|
|
|
|
|
# Landsat too blue-ish because panchromatic band less sensitive to blue
|
|
|
# light, so output blue channed can be modified
|
|
|
if bladjust:
|
|
|
grass.message(_("Adjusting blue channel color table..."))
|
|
|
- blue_colors = ['0 0 0 0\n5% 0 0 0\n67% 255 255 255\n100% 255 255 255']
|
|
|
+ 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(grass.encode('\n'.join(blue_colors)))
|
|
|
+ bc = grass.feed_command("r.colors", quiet=True, map="%s_blue" % out, rules="-")
|
|
|
+ bc.stdin.write(grass.encode("\n".join(blue_colors)))
|
|
|
bc.stdin.close()
|
|
|
|
|
|
# output notice
|
|
|
grass.verbose(_("The following pan-sharpened output maps have been generated:"))
|
|
|
- for ch in ['red', 'green', 'blue']:
|
|
|
+ for ch in ["red", "green", "blue"]:
|
|
|
grass.verbose(_("%s_%s") % (out, ch))
|
|
|
|
|
|
grass.verbose(_("To visualize output, run: g.region -p raster=%s_red" % out))
|
|
|
grass.verbose(_("d.rgb r=%s_red g=%s_green b=%s_blue" % (out, out, out)))
|
|
|
- grass.verbose(_("If desired, combine channels into a single RGB map with 'r.composite'."))
|
|
|
+ grass.verbose(
|
|
|
+ _("If desired, combine channels into a single RGB map with 'r.composite'.")
|
|
|
+ )
|
|
|
grass.verbose(_("Channel colors can be rebalanced using i.colors.enhance."))
|
|
|
|
|
|
# write cmd history:
|
|
|
- for ch in ['red', 'green', 'blue']:
|
|
|
+ for ch in ["red", "green", "blue"]:
|
|
|
grass.raster_history("%s_%s" % (out, ch))
|
|
|
|
|
|
# create a group with the three outputs
|
|
|
- #grass.run_command('i.group', group=out,
|
|
|
+ # grass.run_command('i.group', group=out,
|
|
|
# input="{n}_red,{n}_blue,{n}_green".format(n=out))
|
|
|
|
|
|
# Cleanup
|
|
|
grass.message(_("cleaning up temp files"))
|
|
|
try:
|
|
|
- grass.run_command('g.remove', flags="f", type="raster",
|
|
|
- pattern="tmp%s*" % pid, quiet=True)
|
|
|
+ grass.run_command(
|
|
|
+ "g.remove", flags="f", type="raster", pattern="tmp%s*" % pid, quiet=True
|
|
|
+ )
|
|
|
except:
|
|
|
""
|
|
|
|
|
|
+
|
|
|
def brovey(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
grass.verbose(_("Using Brovey algorithm"))
|
|
|
|
|
|
# pan/intensity histogram matching using linear regression
|
|
|
grass.message(_("Pan channel/intensity histogram matching using linear regression"))
|
|
|
- outname = 'tmp%s_pan1' % pid
|
|
|
+ outname = "tmp%s_pan1" % pid
|
|
|
panmatch1 = matchhist(pan, ms1, outname)
|
|
|
|
|
|
- outname = 'tmp%s_pan2' % pid
|
|
|
+ outname = "tmp%s_pan2" % pid
|
|
|
panmatch2 = matchhist(pan, ms2, outname)
|
|
|
|
|
|
- outname = 'tmp%s_pan3' % pid
|
|
|
+ outname = "tmp%s_pan3" % pid
|
|
|
panmatch3 = matchhist(pan, ms3, outname)
|
|
|
|
|
|
- outr = '%s_red' % out
|
|
|
- outg = '%s_green' % out
|
|
|
- outb = '%s_blue' % out
|
|
|
+ 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")
|
|
|
+ e = """eval(k = "$ms1" + "$ms2" + "$ms3")
|
|
|
"$outr" = 1 * round("$ms3" * "$panmatch3" / k)
|
|
|
"$outg" = 1 * round("$ms2" * "$panmatch2" / k)
|
|
|
- "$outb" = 1 * round("$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)
|
|
|
+ "$outb" = 1 * round("$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 * round((%s * %s) / (%s + %s + %s))' %
|
|
|
- (out, ms1, panmatch1, ms1, ms2, ms3),
|
|
|
- overwrite=True)
|
|
|
- pg = grass.mapcalc_start('%s_green = 1 * round((%s * %s) / (%s + %s + %s))' %
|
|
|
- (out, ms2, panmatch2, ms1, ms2, ms3),
|
|
|
- overwrite=True)
|
|
|
- pr = grass.mapcalc_start('%s_red = 1 * round((%s * %s) / (%s + %s + %s))' %
|
|
|
- (out, ms3, panmatch3, ms1, ms2, ms3),
|
|
|
- overwrite=True)
|
|
|
+ pb = grass.mapcalc_start(
|
|
|
+ "%s_blue = 1 * round((%s * %s) / (%s + %s + %s))"
|
|
|
+ % (out, ms1, panmatch1, ms1, ms2, ms3),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
+ pg = grass.mapcalc_start(
|
|
|
+ "%s_green = 1 * round((%s * %s) / (%s + %s + %s))"
|
|
|
+ % (out, ms2, panmatch2, ms1, ms2, ms3),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
+ pr = grass.mapcalc_start(
|
|
|
+ "%s_red = 1 * round((%s * %s) / (%s + %s + %s))"
|
|
|
+ % (out, ms3, panmatch3, ms1, ms2, ms3),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
|
|
|
pb.wait(), pg.wait(), pr.wait()
|
|
|
try:
|
|
@@ -357,22 +527,31 @@ def brovey(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
|
|
|
# Cleanup
|
|
|
try:
|
|
|
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|
|
- name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|
|
+ grass.run_command(
|
|
|
+ "g.remove",
|
|
|
+ flags="f",
|
|
|
+ quiet=True,
|
|
|
+ type="raster",
|
|
|
+ name="%s,%s,%s" % (panmatch1, panmatch2, panmatch3),
|
|
|
+ )
|
|
|
except:
|
|
|
""
|
|
|
|
|
|
+
|
|
|
def ihs(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
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)
|
|
|
+ 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,
|
|
|
+ )
|
|
|
|
|
|
# pan/intensity histogram matching using linear regression
|
|
|
target = "tmp%s_int" % pid
|
|
@@ -381,45 +560,54 @@ def ihs(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
|
|
|
# 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)
|
|
|
+ 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
|
|
|
try:
|
|
|
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|
|
- name=panmatch)
|
|
|
+ grass.run_command(
|
|
|
+ "g.remove", flags="f", quiet=True, type="raster", name=panmatch
|
|
|
+ )
|
|
|
except:
|
|
|
""
|
|
|
|
|
|
+
|
|
|
def pca(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
|
|
|
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)
|
|
|
+ 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]))
|
|
|
+ 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
|
|
|
+ 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]
|
|
@@ -431,97 +619,87 @@ def pca(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
b3evect3 = b3evect[2]
|
|
|
|
|
|
# Histogram matching
|
|
|
- outname = 'tmp%s_pan1' % pid
|
|
|
+ outname = "tmp%s_pan1" % pid
|
|
|
panmatch1 = matchhist(pan, ms1, outname)
|
|
|
|
|
|
- outname = 'tmp%s_pan2' % pid
|
|
|
+ outname = "tmp%s_pan2" % pid
|
|
|
panmatch2 = matchhist(pan, ms2, outname)
|
|
|
|
|
|
- outname = 'tmp%s_pan3' % pid
|
|
|
+ outname = "tmp%s_pan3" % pid
|
|
|
panmatch3 = matchhist(pan, ms3, outname)
|
|
|
|
|
|
grass.message(_("Performing inverse PCA ..."))
|
|
|
|
|
|
# Get mean value of each channel
|
|
|
- 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'])
|
|
|
+ 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"])
|
|
|
|
|
|
if sproc:
|
|
|
# serial processing
|
|
|
- outr = '%s_red' % out
|
|
|
- outg = '%s_green' % out
|
|
|
- outb = '%s_blue' % out
|
|
|
+ outr = "%s_red" % out
|
|
|
+ outg = "%s_green" % out
|
|
|
+ outb = "%s_blue" % out
|
|
|
|
|
|
cmd1 = "$outb = 1 * round(($panmatch1 * $b1evect1) + ($pca2 * $b1evect2) + ($pca3 * $b1evect3) + $b1mean)"
|
|
|
cmd2 = "$outg = 1 * round(($panmatch2 * $b2evect1) + ($pca2 * $b2evect2) + ($pca3 * $b2evect3) + $b2mean)"
|
|
|
cmd3 = "$outr = 1 * round(($panmatch3 * $b3evect1) + ($pca2 * $b3evect2) + ($pca3 * $b3evect3) + $b3mean)"
|
|
|
|
|
|
- cmd = '\n'.join([cmd1, cmd2, cmd3])
|
|
|
-
|
|
|
- grass.mapcalc(cmd, outb=outb, outg=outg, outr=outr,
|
|
|
- panmatch1=panmatch1,
|
|
|
- panmatch2=panmatch2,
|
|
|
- panmatch3=panmatch3,
|
|
|
- pca2=pca2,
|
|
|
- pca3=pca3,
|
|
|
- b1evect1=b1evect1,
|
|
|
- b2evect1=b2evect1,
|
|
|
- b3evect1=b3evect1,
|
|
|
- b1evect2=b1evect2,
|
|
|
- b2evect2=b2evect2,
|
|
|
- b3evect2=b3evect2,
|
|
|
- b1evect3=b1evect3,
|
|
|
- b2evect3=b2evect3,
|
|
|
- b3evect3=b3evect3,
|
|
|
- b1mean=b1mean,
|
|
|
- b2mean=b2mean,
|
|
|
- b3mean=b3mean,
|
|
|
- overwrite=True)
|
|
|
+ cmd = "\n".join([cmd1, cmd2, cmd3])
|
|
|
+
|
|
|
+ grass.mapcalc(
|
|
|
+ cmd,
|
|
|
+ outb=outb,
|
|
|
+ outg=outg,
|
|
|
+ outr=outr,
|
|
|
+ panmatch1=panmatch1,
|
|
|
+ panmatch2=panmatch2,
|
|
|
+ panmatch3=panmatch3,
|
|
|
+ pca2=pca2,
|
|
|
+ pca3=pca3,
|
|
|
+ b1evect1=b1evect1,
|
|
|
+ b2evect1=b2evect1,
|
|
|
+ b3evect1=b3evect1,
|
|
|
+ b1evect2=b1evect2,
|
|
|
+ b2evect2=b2evect2,
|
|
|
+ b3evect2=b3evect2,
|
|
|
+ b1evect3=b1evect3,
|
|
|
+ b2evect3=b2evect3,
|
|
|
+ b3evect3=b3evect3,
|
|
|
+ b1mean=b1mean,
|
|
|
+ b2mean=b2mean,
|
|
|
+ b3mean=b3mean,
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
else:
|
|
|
# parallel processing
|
|
|
- pb = grass.mapcalc_start('%s_blue = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|
|
- % (out,
|
|
|
- panmatch1,
|
|
|
- b1evect1,
|
|
|
- pca2,
|
|
|
- b1evect2,
|
|
|
- pca3,
|
|
|
- b1evect3,
|
|
|
- b1mean),
|
|
|
- overwrite=True)
|
|
|
-
|
|
|
- pg = grass.mapcalc_start('%s_green = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|
|
- % (out,
|
|
|
- panmatch2,
|
|
|
- b2evect1,
|
|
|
- pca2,
|
|
|
- b2evect2,
|
|
|
- pca3,
|
|
|
- b2evect3,
|
|
|
- b2mean),
|
|
|
- overwrite=True)
|
|
|
-
|
|
|
- pr = grass.mapcalc_start('%s_red = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)'
|
|
|
- % (out,
|
|
|
- panmatch3,
|
|
|
- b3evect1,
|
|
|
- pca2,
|
|
|
- b3evect2,
|
|
|
- pca3,
|
|
|
- b3evect3,
|
|
|
- b3mean),
|
|
|
- overwrite=True)
|
|
|
+ pb = grass.mapcalc_start(
|
|
|
+ "%s_blue = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)"
|
|
|
+ % (out, panmatch1, b1evect1, pca2, b1evect2, pca3, b1evect3, b1mean),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
+
|
|
|
+ pg = grass.mapcalc_start(
|
|
|
+ "%s_green = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)"
|
|
|
+ % (out, panmatch2, b2evect1, pca2, b2evect2, pca3, b2evect3, b2mean),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
+
|
|
|
+ pr = grass.mapcalc_start(
|
|
|
+ "%s_red = 1 * round((%s * %f) + (%s * %f) + (%s * %f) + %f)"
|
|
|
+ % (out, panmatch3, b3evect1, pca2, b3evect2, pca3, b3evect3, b3mean),
|
|
|
+ overwrite=True,
|
|
|
+ )
|
|
|
|
|
|
pb.wait(), pg.wait(), pr.wait()
|
|
|
try:
|
|
@@ -530,16 +708,22 @@ def pca(pan, ms1, ms2, ms3, out, pid, sproc):
|
|
|
""
|
|
|
|
|
|
# Cleanup
|
|
|
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|
|
- name='%s,%s,%s' % (panmatch1, panmatch2, panmatch3))
|
|
|
+ grass.run_command(
|
|
|
+ "g.remove",
|
|
|
+ flags="f",
|
|
|
+ quiet=True,
|
|
|
+ type="raster",
|
|
|
+ name="%s,%s,%s" % (panmatch1, panmatch2, panmatch3),
|
|
|
+ )
|
|
|
+
|
|
|
|
|
|
def matchhist(original, target, matched):
|
|
|
# pan/intensity histogram matching using numpy arrays
|
|
|
grass.message(_("Histogram matching..."))
|
|
|
|
|
|
# input images
|
|
|
- original = original.split('@')[0]
|
|
|
- target = target.split('@')[0]
|
|
|
+ original = original.split("@")[0]
|
|
|
+ target = target.split("@")[0]
|
|
|
images = [original, target]
|
|
|
|
|
|
# create a dictionary to hold arrays for each image
|
|
@@ -547,14 +731,13 @@ def matchhist(original, target, matched):
|
|
|
|
|
|
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=img,
|
|
|
- sep=':')
|
|
|
- stats = grass.decode(stats_out.communicate()[0]).split('\n')[:-1]
|
|
|
- stats_dict = dict(s.split(':', 1) for s in stats)
|
|
|
+ stats_out = grass.pipe_command("r.stats", flags="cin", input=img, sep=":")
|
|
|
+ stats = grass.decode(stats_out.communicate()[0]).split("\n")[:-1]
|
|
|
+ stats_dict = dict(s.split(":", 1) for s in stats)
|
|
|
total_cells = 0 # total non-null cells
|
|
|
for j in stats_dict:
|
|
|
stats_dict[j] = int(stats_dict[j])
|
|
|
- if j != '*':
|
|
|
+ if j != "*":
|
|
|
total_cells += stats_dict[j]
|
|
|
|
|
|
if total_cells < 1:
|
|
@@ -564,8 +747,10 @@ def matchhist(original, target, matched):
|
|
|
# cumulative distribution function (CDF) for each grey value.
|
|
|
# Grey value is the integer (i4) and cdf is float (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
|
|
|
+ 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):
|
|
|
if str(n) in stats_dict:
|
|
@@ -583,7 +768,7 @@ def matchhist(original, target, matched):
|
|
|
arrays[img][n] = (n, cdf)
|
|
|
|
|
|
# open file for reclass rules
|
|
|
- outfile = open(grass.tempfile(), 'w')
|
|
|
+ outfile = open(grass.tempfile(), "w")
|
|
|
|
|
|
for i in arrays[original]:
|
|
|
# for each grey value and corresponding cdf value in original, find the
|
|
@@ -610,15 +795,14 @@ def matchhist(original, target, matched):
|
|
|
outfile.close()
|
|
|
|
|
|
# create reclass of target from reclass rules file
|
|
|
- result = grass.core.find_file(matched, element='cell')
|
|
|
- if result['fullname']:
|
|
|
- grass.run_command('g.remove', flags='f', quiet=True, type='raster',
|
|
|
- name=matched)
|
|
|
- grass.run_command('r.reclass', input=original, out=matched,
|
|
|
- rules=outfile.name)
|
|
|
+ result = grass.core.find_file(matched, element="cell")
|
|
|
+ if result["fullname"]:
|
|
|
+ grass.run_command(
|
|
|
+ "g.remove", flags="f", quiet=True, type="raster", name=matched
|
|
|
+ )
|
|
|
+ grass.run_command("r.reclass", input=original, out=matched, rules=outfile.name)
|
|
|
else:
|
|
|
- grass.run_command('r.reclass', input=original, out=matched,
|
|
|
- rules=outfile.name)
|
|
|
+ grass.run_command("r.reclass", input=original, out=matched, rules=outfile.name)
|
|
|
|
|
|
# Cleanup
|
|
|
# remove the rules file
|
|
@@ -627,6 +811,7 @@ def matchhist(original, target, matched):
|
|
|
# return reclass of target with histogram that matches original
|
|
|
return matched
|
|
|
|
|
|
+
|
|
|
if __name__ == "__main__":
|
|
|
options, flags = grass.parser()
|
|
|
main()
|