Convert d.correlate to Python

Add, use feed_command, write_command (write to stdin)
Clean up r.in.aster


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@33562 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/python/grass.py

@@ -54,10 +54,22 @@ def pipe_command(*args, **kwargs):
     kwargs['stdout'] = subprocess.PIPE
     return start_command(*args, **kwargs)
 
+def feed_command(*args, **kwargs):
+    kwargs['stdin'] = subprocess.PIPE
+    return start_command(*args, **kwargs)
+
 def read_command(*args, **kwargs):
     ps = pipe_command(*args, **kwargs)
     return ps.communicate()[0]
 
+def write_command(*args, **kwargs):
+    stdin = kwargs['stdin']
+    kwargs['stdin'] = subprocess.PIPE
+    p = start_command(*args, **kwargs)
+    p.stdin.write(stdin)
+    p.stdin.close()
+    return p.wait()
+
 def exec_command(prog, flags = "", overwrite = False, quiet = False, verbose = False, env = None, **kwargs):
     args = make_command(prog, flags, overwrite, quiet, verbose, **kwargs)
     if env == None:
@@ -274,10 +286,10 @@ def basename(path, ext = None):
 
 # find a program (replacement for "which")
 
-def find_program(pgm):
+def find_program(pgm, args = []):
     nuldev = file(os.devnull, 'w+')
     try:
-	subprocess.call([pgm], stdin = nuldev, stdout = nuldev, stderr = nuldev)
+	subprocess.call([pgm] + args, stdin = nuldev, stdout = nuldev, stderr = nuldev)
 	found = True
     except:
 	found = False

scripts/d.correlate/d.correlate.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+#
+############################################################################
+#
+# MODULE:       d.correlate for GRASS 6; based on dcorrelate.sh for GRASS 4,5
+# AUTHOR(S):    CERL - Michael Shapiro; updated to GRASS 6 by Markus Neteler 5/2005
+#               Converted to Python by Glynn Clements
+# PURPOSE:      prints a graph of the correlation between data layers (in pairs)
+#               derived from <grass5>/src.local/d.correlate.sh
+# COPYRIGHT:    (C) 2005, 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.
+#
+# TODO GRASS 7: rename parameters to map1, map2 or better multiple map=map1[,map2[...]]
+#############################################################################
+
+#%Module
+#%  description: Prints a graph of the correlation between data layers (in pairs).
+#%  keywords: display, diagram
+#%End
+#%option
+#% key: layers
+#% type: string
+#% gisprompt: old,cell,raster
+#% description: raster input map
+#% required: yes
+#% multiple: yes
+#%end
+
+import sys
+import os
+import grass
+
+def main():
+    layers = options['layers'].split(',')
+
+    if len(layers) < 2:
+	grass.error("At least 2 layers are required")
+
+    tmpfile = grass.tempfile()
+
+    for map in layers:
+	if not grass.find_file(map, element = 'cell')['file']:
+	    grass.fatal("Input map <%s> not found" % map)
+
+    grass.write_command('d.text', color = 'black', size = 4, line = 1, stdin = "CORRELATION")
+
+    colors = "red black blue green gray violet".split()
+    line = 2
+    iloop = 0
+    jloop = 0
+    for iloop, i in enumerate(layers):
+	for jloop, j in enumerate(layers):
+	    if i != j and iloop <= jloop:
+		color = colors[0]
+		colors = colors[1:]
+		colors.append(color)
+		grass.write_command('d.text', color = color, size = 4, line = line, stdin = "%s %s" % (i, j))
+		line += 1
+
+		ofile = file(tmpfile, 'w')
+		grass.run_command('r.stats', flags = 'cnA', input = (i, j), stdout = ofile)
+		ofile.close()
+
+		ifile = file(tmpfile, 'r')
+		first = True
+		for l in ifile:
+		    f = l.rstrip('\r\n').split(' ')
+		    x = float(f[0])
+		    y = float(f[1])
+		    if first:
+			minx = maxx = x
+			miny = maxy = y
+			first = False
+		    if minx > x: minx = x
+		    if maxx < x: maxx = x
+		    if miny > y: miny = y
+		    if maxy < y: maxy = y
+		ifile.close()
+
+		kx = 100.0/(maxx-minx+1)
+		ky = 100.0/(maxy-miny+1)
+
+		p = grass.feed_command('d.graph', color = color)
+		ofile = p.stdin
+
+		ifile = file(tmpfile, 'r')
+		for l in ifile:
+		    f = l.rstrip('\r\n').split(' ')
+		    x = float(f[0])
+		    y = float(f[1])
+		    ofile.write("icon + 0.1 %f %f\n" % ((x-minx+1) * kx, (y-miny+1) * ky))
+		ifile.close()
+
+		ofile.close()
+		p.wait()
+
+    grass.try_remove(tmpfile)
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()

scripts/db.droptable/db.droptable.py

@@ -37,7 +37,6 @@
 import sys
 import os
 import grass
-import subprocess
 
 def main():
     table = options['table']
@@ -84,8 +83,7 @@ def main():
 	grass.message("You must use the force flag to actually remove it. Exiting.")
 	sys.exit(0)
 
-    p = grass.start_command('db.execute', database = database, driver = driver,
-			    stdin = subprocess.PIPE)
+    p = grass.feed_command('db.execute', database = database, driver = driver)
     p.stdin.write("DROP TABLE " + table)
     p.stdin.close()
     p.wait()

scripts/r.in.aster/r.in.aster.py

@@ -101,37 +101,15 @@ bands = {
     }
 }
 
-def _message(msg, *args):
-    subprocess.call(["g.message", "message=%s" % msg] + list(args))
-
-def debug(msg):
-    _message(msg, '-d')
-
-def message(msg):
-    _message(msg)
-
-def error(msg):
-    _message(msg, '-e')
-    sys.exit(1)
-
 def main():
+    input = options['input']
+    proctype = options['proctype']
+    output = options['output']
+    band = options['band']
 
     #check whether gdalwarp is in path and executable
-    p = None
-    try:
-        p = subprocess.call(['gdalwarp', '--version'])
-    except:
-        pass
-    if p == None or p != 0:
-        error("gdalwarp is not in the path and executable")
-
-    #initialize variables
-    dataset = ''
-    srcfile = ''
-    proj = ''
-    band = ''
-    outfile = ''
-    bandlist = []
+    if not grass.find_command('gdalwarp', ['--version']):
+        grass.fatal("gdalwarp is not in the path and executable")
 
     #create temporary file to hold gdalwarp output before importing to GRASS
     tempfile = grass.read_command("g.tempfile", pid = os.getpid()).strip() + '.tif'
@@ -141,41 +119,41 @@ def main():
 
     #currently only runs in projected location
     if "XY location" in proj:
-      error ("This module needs to be run in a projected location (found: %s)" % proj)
+      grass.fatal("This module needs to be run in a projected location (found: %s)" % proj)
 
 
     #process list of bands
     allbands = ['1','2','3n','3b','4','5','6','7','8','9','10','11','12','13','14']
-    if options['band'].strip() == 'all':
+    if band == 'all':
         bandlist = allbands
     else:
-        bandlist = options['band'].split(',')
+        bandlist = band.split(',')
 
     #initialize datasets for L1A and L1B
-    if options['proctype'] in ["L1A", "L1B"]:
+    if proctype in ["L1A", "L1B"]:
         for band in bandlist:
             if band in allbands:
-                dataset = bands[options['proctype']][band]
-                srcfile = "HDF4_EOS:EOS_SWATH:%s:%s" % (options['input'], dataset)
+                dataset = bands[proctype][band]
+                srcfile = "HDF4_EOS:EOS_SWATH:%s:%s" % (input, dataset)
                 import_aster(proj, srcfile, tempfile, band)
             else:
-                error('band %s is not an available Terra/ASTER band' % band)
-    elif options['proctype'] == "DEM": 
-        srcfile=options['input']
+                grass.fatal('band %s is not an available Terra/ASTER band' % band)
+    elif proctype == "DEM": 
+        srcfile = input
         import_aster(proj, srcfile, tempfile, "DEM")
 
     #cleanup
-    message("Cleaning up ...")
+    grass.message("Cleaning up ...")
     grass.try_remove(tempfile)
-    message("Done.")
+    grass.message("Done.")
 
     return
 
 def import_aster(proj, srcfile, tempfile, band):
     #run gdalwarp with selected options (must be in $PATH)
     #to translate aster image to geotiff
-    message("Georeferencing aster image ...")
-    debug("gdalwarp -t_srs %s %s %s" % (proj, srcfile, tempfile))
+    grass.message("Georeferencing aster image ...")
+    grass.debug("gdalwarp -t_srs %s %s %s" % (proj, srcfile, tempfile))
 
     if platform.system() == "Darwin":
         cmd = ["arch", "-i386", "gdalwarp", "-t_srs", proj, srcfile, tempfile ]
@@ -189,8 +167,8 @@ def import_aster(proj, srcfile, tempfile, band):
     #p = subprocess.call(["gdal_translate", srcfile, tempfile])
 
     #import geotiff to GRASS
-    message("Importing into GRASS ...")
-    outfile = options['output'].strip()+'.'+band
+    grass.message("Importing into GRASS ...")
+    outfile = "%s.%s" % (output, band)
     grass.run_command("r.in.gdal", overwrite = flags['o'], input = tempfile, output = outfile)
 
     # write cmd history

scripts/r.mask/r.mask.py

@@ -45,7 +45,6 @@ import sys
 import os
 import grass
 import atexit
-import subprocess
 
 def cleanup():
     if tmp:
@@ -73,8 +72,7 @@ def main():
 	if exists and not grass.overwrite():
 	    grass.fatal("MASK already found in current mapset. Delete first or overwrite")
 
-	p = grass.start_command('r.reclass', input = input, output = 'MASK', overwrite = True,
-				stdin = subprocess.PIPE)
+	p = grass.feed_command('r.reclass', input = input, output = 'MASK', overwrite = True)
 	p.stdin.write("%s = 1" % maskcats)
 	p.stdin.close()
 	p.wait()

scripts/v.db.addcol/v.db.addcol.py

@@ -48,7 +48,6 @@
 import sys
 import os
 import grass
-import subprocess
 
 def main():
     map = options['map']
@@ -77,7 +76,7 @@ def main():
 	if not col:
 	    grass.fatal("There is an empty column. Did you leave a trailing comma?")
 
-	p = grass.start_command('db.execute', database = database, driver = driver, stdin = subprocess.PIPE)
+	p = grass.feed_command('db.execute', database = database, driver = driver)
 	p.stdin.write("ALTER TABLE %s ADD COLUMN %s" % (table, col))
 	p.stdin.close()
 	if p.wait() != 0:

scripts/v.db.join/v.db.join.py

@@ -59,7 +59,6 @@
 
 import sys
 import os
-import subprocess
 import string
 import grass
 
@@ -112,7 +111,7 @@ def main():
 				   otable = otable, ocolumn = ocolumn,
 				   colname = colname)
 
-	p = grass.start_command('db.execute', stdin = subprocess.PIPE)
+	p = grass.feed_command('db.execute')
 	p.stdin.write(stmt)
 	p.stdin.close()
 	p.wait()