[GRASS-SVN] r33632 - in grass/trunk/scripts: i.fusion.brovey
i.image.mosaic i.landsat.rgb r.in.aster r.shaded.relief
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Sep 30 20:58:59 EDT 2008
Author: glynn
Date: 2008-09-30 20:58:59 -0400 (Tue, 30 Sep 2008)
New Revision: 33632
Added:
grass/trunk/scripts/i.fusion.brovey/i.fusion.brovey.py
Modified:
grass/trunk/scripts/i.image.mosaic/i.image.mosaic.py
grass/trunk/scripts/i.landsat.rgb/i.landsat.rgb.py
grass/trunk/scripts/r.in.aster/r.in.aster.py
grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
Log:
Convert i.fusion.brovey to Python
Fix various bugs in other scripts
Added: grass/trunk/scripts/i.fusion.brovey/i.fusion.brovey.py
===================================================================
--- grass/trunk/scripts/i.fusion.brovey/i.fusion.brovey.py (rev 0)
+++ grass/trunk/scripts/i.fusion.brovey/i.fusion.brovey.py 2008-10-01 00:58:59 UTC (rev 33632)
@@ -0,0 +1,185 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: i.fusion.brovey
+# AUTHOR(S): Markus Neteler. <neteler itc it>
+# Converted to Python by Glynn Clements
+# PURPOSE: Brovey transform to merge
+# - LANDSAT-7 MS (2, 4, 5) and pan (high res)
+# - SPOT MS and pan (high res)
+# - QuickBird MS and pan (high res)
+#
+# COPYRIGHT: (C) 2002-2008 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.
+#
+# REFERENCES:
+# (?) Roller, N.E.G. and Cox, S., 1980. Comparison of Landsat MSS
+# and merged MSS/RBV data for analysis of natural vegetation.
+# Proc. of the 14th International Symposium on Remote Sensing
+# of Environment, San Jose, Costa Rica, 23-30 April, pp. 1001-1007.
+#
+# for LANDSAT 5: see Pohl, C 1996 and others
+#
+# TODO: add overwrite test at beginning of the script
+#############################################################################
+
+#%Module
+#% description: Brovey transform to merge multispectral and high-res panchromatic channels
+#% keywords: raster, imagery, fusion
+#%End
+#%Flag
+#% key: l
+#% description: sensor: LANDSAT
+#%END
+#%Flag
+#% key: q
+#% description: sensor: QuickBird
+#%END
+#%Flag
+#% key: s
+#% description: sensor: SPOT
+#%END
+#%option
+#% key: ms1
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: raster input map (green: tm2 | qbird_green | spot1)
+#% required : yes
+#%end
+#%option
+#% key: ms2
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: raster input map (NIR: tm4 | qbird_nir | spot2)
+#% required : yes
+#%end
+#%option
+#% key: ms3
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: raster input map (MIR; tm5 | qbird_red | spot3)
+#% required : yes
+#%end
+#%option
+#% key: pan
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: raster input map (etmpan | qbird_pan | spotpan)
+#% required : yes
+#%end
+#%option
+#% key: outputprefix
+#% type: string
+#% gisprompt: new,cell,raster
+#% description: raster output map prefix (e.g. 'brov')
+#% required : yes
+#%end
+
+import sys
+import os
+import string
+import grass
+
+def main():
+ global tmp
+
+ landsat = flags['l']
+ quickbird = flags['q']
+ spot = flags['s']
+
+ ms1 = options['ms1']
+ ms2 = options['ms2']
+ ms3 = options['ms3']
+ pan = options['pan']
+ out = options['outputprefix']
+
+ tmp = str(os.getpid())
+
+ if not landsat and not quickbird and not spot:
+ grass.fatal("Please select a flag to specify the satellite sensor")
+
+ #get PAN resolution:
+ s = grass.read_command('r.info', flags = 's', map = pan)
+ kv = grass.parse_key_val(s)
+ nsres = float(kv['nsres'])
+ ewres = float(kv['ewres'])
+ panres = (nsres + ewres) / 2
+
+ # clone current region
+ grass.use_temp_region()
+
+ grass.message("Using resolution from PAN: %f" % panres)
+ grass.run_command('g.region', flags = 'a', res = panres)
+
+ grass.message("Performing Brovey transformation...")
+
+ # The formula was originally developed for LANDSAT-TM5 and SPOT,
+ # but it also works well with LANDSAT-TM7
+ # LANDSAT formula:
+ # r.mapcalc "brov.red=1. * tm.5 / (tm.2 + tm.4 + tm.5) * etmpan"
+ # r.mapcalc "brov.green=1. * tm.4 /(tm.2 + tm.4 + tm.5) * etmpan"
+ # r.mapcalc "brov.blue=1. * tm.2 / (tm.2 + tm.4 + tm.5) * etmpan"
+ #
+ # SPOT formula:
+ # r.mapcalc "brov.red= 1. * spot.ms.3 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
+ # r.mapcalc "brov.green=1. * spot.ms.2 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
+ # r.mapcalc "brov.blue= 1. * spot.ms.1 / (spot.ms.1 + spot.ms.2 + spot.ms.3) * spot.p"
+ # note: for RGB composite then revert brov.red and brov.green!
+
+ grass.message("Calculating %s.{red,green,blue}: ..." % out)
+ t = string.Template(
+ '''"$out.red" = eval(k = "$pan" / ("$ms1" + "$ms2" + "$ms3"), 1. * "$ms3" * k)
+ "$out.green" = 1. * "$ms2" * k
+ "$out.blue" = 1. * "$ms1" * k''')
+ e = t.substitute(out = out, pan = pan, ms1 = ms1, ms2 = ms2, ms3 = ms3)
+ if grass.run_command('r.mapcalc', expression = e) != 0:
+ grass.fatal("An error occurred while running r.mapcalc")
+
+ # Maybe?
+ #r.colors $GIS_OPT_OUTPUTPREFIX.red col=grey
+ #r.colors $GIS_OPT_OUTPUTPREFIX.green col=grey
+ #r.colors $GIS_OPT_OUTPUTPREFIX.blue col=grey
+ #to blue-ish, therefore we modify
+ #r.colors $GIS_OPT_OUTPUTPREFIX.blue col=rules << EOF
+ #5 0 0 0
+ #20 200 200 200
+ #40 230 230 230
+ #67 255 255 255
+ #EOF
+
+ if spot:
+ #apect table is nice for SPOT:
+ grass.message("Assigning color tables for SPOT ...")
+ for ch in ['red', 'green', 'blue']:
+ grass.run_command('r.colors', map = "%s.%s" % (out, ch), col = 'aspect')
+ grass.message("Fixing output names...")
+ for s, d in [('green','tmp'),('red','green'),('tmp','red')]:
+ src = "%s.%s" % (out, s)
+ dst = "%s.%s" % (out, d)
+ grass.run_command('g.rename', rast = (src, dst), quiet = True)
+ else:
+ #aspect table is nice for LANDSAT and QuickBird:
+ grass.message("Assigning color tables for LANDSAT or QuickBird ...")
+ for ch in ['red', 'green', 'blue']:
+ grass.run_command('r.colors', map = "%s.%s" % (out, ch), col = 'aspect')
+
+ grass.message("Following pan-sharpened output maps have been generated:")
+ for ch in ['red', 'green', 'blue']:
+ grass.message("%s.%s" % (out, ch))
+
+ grass.message("To visualize output, run:")
+ grass.message("g.region -p rast=%s.red" % out)
+ grass.message("d.rgb r=%s.red g=%s.green b=%s.blue" % (out, out, out))
+ grass.message("If desired, combine channels with 'r.composite' to a single map.")
+
+ # write cmd history:
+ for ch in ['red', 'green', 'blue']:
+ grass.raster_history("%s.%s" % (out, ch))
+
+if __name__ == "__main__":
+ options, flags = grass.parser()
+ main()
Property changes on: grass/trunk/scripts/i.fusion.brovey/i.fusion.brovey.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: grass/trunk/scripts/i.image.mosaic/i.image.mosaic.py
===================================================================
--- grass/trunk/scripts/i.image.mosaic/i.image.mosaic.py 2008-10-01 00:57:14 UTC (rev 33631)
+++ grass/trunk/scripts/i.image.mosaic/i.image.mosaic.py 2008-10-01 00:58:59 UTC (rev 33632)
@@ -108,7 +108,7 @@
grass.message("Ready. File %s created." % output)
# write cmd history:
- raster_history(output)
+ grass.raster_history(output)
if __name__ == "__main__":
options, flags = grass.parser()
Modified: grass/trunk/scripts/i.landsat.rgb/i.landsat.rgb.py
===================================================================
--- grass/trunk/scripts/i.landsat.rgb/i.landsat.rgb.py 2008-10-01 00:57:14 UTC (rev 33631)
+++ grass/trunk/scripts/i.landsat.rgb/i.landsat.rgb.py 2008-10-01 00:58:59 UTC (rev 33632)
@@ -132,7 +132,7 @@
# write cmd history:
for i in [red, green, blue]:
- raster_history(i)
+ grass.raster_history(i)
if __name__ == "__main__":
options, flags = grass.parser()
Modified: grass/trunk/scripts/r.in.aster/r.in.aster.py
===================================================================
--- grass/trunk/scripts/r.in.aster/r.in.aster.py 2008-10-01 00:57:14 UTC (rev 33631)
+++ grass/trunk/scripts/r.in.aster/r.in.aster.py 2008-10-01 00:58:59 UTC (rev 33632)
@@ -172,7 +172,7 @@
grass.run_command("r.in.gdal", overwrite = flags['o'], input = tempfile, output = outfile)
# write cmd history
- raster_history(outfile)
+ grass.raster_history(outfile)
if __name__ == "__main__":
options, flags = grass.parser()
Modified: grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py
===================================================================
--- grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2008-10-01 00:57:14 UTC (rev 33631)
+++ grass/trunk/scripts/r.shaded.relief/r.shaded.relief.py 2008-10-01 00:58:59 UTC (rev 33632)
@@ -117,8 +117,9 @@
azimuth = options['azimuth']
zmult = options['zmult']
scale = options['scale']
+ units = options['units']
- if not grass.findfile(elev)['file']:
+ if not grass.find_file(input)['file']:
grass.fatal("Map <%s> not found." % elev)
if input == output:
@@ -147,22 +148,22 @@
scale *= 6076.12 * 60
#correct azimuth to East (GRASS convention):
- az = azimuth - 90
+ az = float(azimuth) - 90
grass.message("Calculating shading, please stand by.")
- t = string.Template(r'$elev_out = eval(
+ t = string.Template(r'''$elev_out = eval( \
x=($zmult*$elev[-1,-1] + 2*$zmult*$elev[0,-1] + $zmult*$elev[1,-1] - \
- $zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale),
+ $zmult*$elev[-1,1] - 2*$zmult*$elev[0,1] - $zmult*$elev[1,1])/(8.*ewres()*$scale), \
y=($zmult*$elev[-1,-1] + 2*$zmult*$elev[-1,0] + $zmult*$elev[-1,1] - \
- $zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale),
-slope=90.-atan(sqrt(x*x + y*y)),
-a=round(atan(x,y)),
-a=if(isnull(a),1,a),
-aspect=if(x!=0||y!=0,if(a,a,360.)),
-cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect),
-if(cang < 0.,0.,100.*cang),
-if(isnull(cang), null(), 100.*cang))')
+ $zmult*$elev[1,-1] - 2*$zmult*$elev[1,0] - $zmult*$elev[1,1])/(8.*nsres()*$scale), \
+slope=90.-atan(sqrt(x*x + y*y)), \
+a=round(atan(x,y)), \
+a=if(isnull(a),1,a), \
+aspect=if(x!=0||y!=0,if(a,a,360.)), \
+cang = sin($alt)*sin(slope) + cos($alt)*cos(slope) * cos($az-aspect), \
+if(cang < 0.,0.,100.*cang), \
+if(isnull(cang), null(), 100.*cang))''')
expr = t.substitute(alt = alt, az = az, elev = elev, elev_out = elev_out, scale = scale, zmult = zmult)
p = grass.feed_command('r.mapcalc')
p.stdin.write(expr)
@@ -188,7 +189,7 @@
grass.message("Shaded relief map created and named <%s>. Consider renaming." % elevshade)
# write cmd history:
- raster_history(elevshade)
+ grass.raster_history(elevshade)
if __name__ == "__main__":
options, flags = grass.parser()
More information about the grass-commit
mailing list