Add the new GridModule class

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

lib/python/pygrass/modules/Makefile

@@ -11,7 +11,7 @@ DSTDIR= $(PGDIR)/modules
 
 MODULES = shortcuts
 
-CLEAN_SUBDIRS = interface
+CLEAN_SUBDIRS = interface grid
 
 PYFILES := $(patsubst %,$(DSTDIR)/%.py,$(MODULES) __init__)
 PYCFILES := $(patsubst %,$(DSTDIR)/%.pyc,$(MODULES) __init__)

lib/python/pygrass/modules/__init__.py

@@ -8,7 +8,7 @@ Created on Thu Jul 12 10:23:15 2012
 import interface
 from interface import Module
 
-#import grid
-#from grid.grid import GridModule
+import grid
+from grid.grid import GridModule
 
 import shortcuts

lib/python/pygrass/modules/grid/Makefile

@@ -0,0 +1,32 @@
+MODULE_TOPDIR = ../../../../..
+
+include $(MODULE_TOPDIR)/include/Make/Other.make
+include $(MODULE_TOPDIR)/include/Make/Python.make
+include $(MODULE_TOPDIR)/include/Make/Doxygen.make
+
+PYDIR = $(ETC)/python
+GDIR = $(PYDIR)/grass
+PGDIR = $(GDIR)/pygrass
+DSTDIR= $(PGDIR)/modules/grid
+
+MODULES = split patch grid
+
+PYFILES := $(patsubst %,$(DSTDIR)/%.py,$(MODULES) __init__)
+PYCFILES := $(patsubst %,$(DSTDIR)/%.pyc,$(MODULES) __init__)
+
+default: $(PYFILES) $(PYCFILES) $(GDIR)/__init__.py $(GDIR)/__init__.pyc
+
+$(PYDIR):
+	$(MKDIR) $@
+
+$(GDIR): | $(PYDIR)
+	$(MKDIR) $@
+
+$(DSTDIR): | $(GDIR)
+	$(MKDIR) $@
+
+$(DSTDIR)/%: % | $(DSTDIR)
+	$(INSTALL_DATA) $< $@
+
+#doxygen:
+DOXNAME = pythonpygrass

lib/python/pygrass/modules/grid/__init__.py

@@ -0,0 +1,10 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr  3 10:41:30 2013
+
+@author: pietro
+"""
+import split
+import patch
+import grid
+from grid import GridModule

lib/python/pygrass/modules/grid/grid.py

@@ -0,0 +1,182 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Mar 28 11:06:00 2013
+
+@author: pietro
+"""
+import os
+import multiprocessing as mltp
+import subprocess as sub
+
+from grass.script.setup import write_gisrc
+
+from grass.pygrass.gis import Mapset, Location, make_mapset
+from grass.pygrass.modules import Module
+
+from split import split_region_tiles
+from patch import patch_map
+
+
+_GREG = Module('g.region')
+
+
+def get_cmd(cmdd):
+    """Transforma a cmd dictionary to a list of parameters"""
+    cmd = [cmdd['name'], ]
+    cmd.extend(("%s=%s" % (k, v) for k, v in cmdd['inputs']))
+    cmd.extend(("%s=%s" % (k, v) for k, v in cmdd['outputs']))
+    cmd.extend(("%s" % (flg) for flg in cmdd['flags'] if len(flg) == 1))
+    cmd.extend(("--%s" % (flg[0]) for flg in cmdd['flags'] if len(flg) > 1))
+    return cmd
+
+
+def cmd_exe((bbox, mapnames, msetname, cmd)):
+    """Create a mapset, and execute a cmd inside."""
+    mset = Mapset()
+    try:
+        make_mapset(msetname)
+    except:
+        pass
+    env = os.environ.copy()
+    env['GISRC'] = write_gisrc(mset.gisdbase, mset.location, msetname)
+    if mapnames:
+        inputs = dict(cmd['inputs'])
+        # reset the inputs to
+        for key in mapnames:
+            inputs[key] = mapnames[key]
+        cmd['inputs'] = inputs.items()
+        # set the region to the tile
+        _GREG(env_=env, rast=key)
+    else:
+        # set the computational region
+        _GREG(env_=env, **bbox)
+    # run the grass command
+    sub.Popen(get_cmd(cmd), env=env).wait()
+
+
+class GridModule(object):
+    """Run GRASS raster commands in a multiproccessing mode.
+
+    Parameters
+    -----------
+
+    cmd: raster GRASS command
+        Only command staring with r.* are valid.
+    width: integer
+        Width of the tile, in pixel.
+    height: integer
+        Height of the tile, in pixel.
+    overlap: integer
+        Overlap between tiles, in pixel.
+    nthreads: number of threads
+        Default value is equal to the number of processor available.
+    split: boolean
+        If True use r.tile to split all the inputs.
+    run_: boolean
+        If False only instantiate the object.
+    args and kargs: cmd parameters
+        Give all the parameters to the command.
+
+    Examples
+    --------
+
+    ::
+
+        >>> grd = GridModule('r.slope.aspect',
+        ...                  width=500, height=500, overlap=2,
+        ...                  processes=None, split=True,
+        ...                  elevation='elevation',
+        ...                  slope='slope', aspect='aspect', overwrite=True)
+        >>> grd.run()
+    """
+    def __init__(self, cmd, width=None, height=None, overlap=0, processes=None,
+                 split=False, *args, **kargs):
+        kargs['run_'] = False
+        self.mset = Mapset()
+        self.module = Module(cmd, *args, **kargs)
+        self.width = width
+        self.height = height
+        self.overlap = overlap
+        self.processes = processes
+        self.bboxes = split_region_tiles(width=width, height=height,
+                                         overlap=overlap)
+        self.msetstr = cmd.replace('.', '') + "_%03d_%03d"
+        self.inlist = None
+        if split:
+            self.split()
+
+    def clean_location(self):
+        """Remove all created mapsets."""
+        mapsets = Location().mapsets(self.msetstr.split('_')[0] + '_*')
+        for mset in mapsets:
+            Mapset(mset).delete()
+
+    def split(self):
+        """Split all the raster inputs using r.tile"""
+        rtile = Module('r.tile')
+        inlist = {}
+        #import pdb; pdb.set_trace()
+        for inmap in self.module.inputs:
+            inm = self.module.inputs[inmap]
+            if inm.type == 'raster' and inm.value:
+                rtile(input=inm.value, output=inm.value,
+                      width=self.width, height=self.height,
+                      overlap=self.overlap)
+                patt = '%s-*' % inm.value
+                inlist[inm.value] = sorted(self.mset.glist(type='rast',
+                                                           pattern=patt))
+        self.inlist = inlist
+
+    def get_works(self):
+        """Return a list of tuble with the parameters for cmd_exe function"""
+        works = []
+        cmd = self.module.get_dict()
+        for row, box_row in enumerate(self.bboxes):
+            for col, box in enumerate(box_row):
+                inms = None
+                if self.inlist:
+                    inms = {}
+                    cols = len(box_row)
+                    for key in self.inlist:
+                        indx = row * cols + col
+                        inms[key] = "%s@%s" % (self.inlist[key][indx],
+                                               self.mset.name)
+                # set the computational region, prepare the region parameters
+                bbox = dict([(k[0], str(v)) for k, v in box.items()[:-2]])
+                works.append((bbox, inms, self.msetstr % (row, col), cmd))
+        return works
+
+    def run(self, patch=True, clean=True):
+        """Run the GRASS command."""
+        self.module.flags.overwrite = True
+        pool = mltp.Pool(processes=self.processes)
+        result = pool.map_async(cmd_exe, self.get_works())
+        result.wait()
+        if not result.successful():
+            raise RuntimeError
+
+        if patch:
+            self.patch()
+
+        if clean:
+            self.clean_location()
+            self.rm_tiles()
+
+    def patch(self):
+        """Patch the final results."""
+        # patch all the outputs
+        for otmap in self.module.outputs:
+            otm = self.module.outputs[otmap]
+            if otm.type == 'raster' and otm.value:
+                patch_map(otm.value, self.mset.name, self.msetstr,
+                          split_region_tiles(width=self.width,
+                                             height=self.height),
+                          self.module.flags.overwrite)
+
+    def rm_tiles(self):
+        """Remove all the tiles."""
+        # if split, remove tiles
+        if self.inlist:
+            grm = Module('g.remove')
+            for key in self.inlist:
+                grm(rast=self.inlist[key])

lib/python/pygrass/modules/grid/patch.py

@@ -0,0 +1,60 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  2 18:57:42 2013
+
+@author: pietro
+"""
+from grass.pygrass.gis.region import Region
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.functions import coor2pixel
+
+
+def get_start_end_index(bbox_list):
+    """Convert a Bounding Box to a list of the index of
+    column start, end, row start and end
+    """
+    ss_list = []
+    reg = Region()
+    for bbox in bbox_list:
+        r_start, c_start = coor2pixel((bbox.west, bbox.north), reg)
+        r_end, c_end = coor2pixel((bbox.east, bbox.south), reg)
+        ss_list.append((int(r_start), int(r_end), int(c_start), int(c_end)))
+    return ss_list
+
+
+def patch_row(rast, rasts, bboxes):
+    """Patch a row of bound boxes."""
+    sei = get_start_end_index(bboxes)
+    # instantiate two buffer
+    buff = rasts[0][0]
+    rbuff = rasts[0][0]
+    r_start, r_end, c_start, c_end = sei[0]
+    for row in xrange(r_start, r_end):
+        for col, ras in enumerate(rasts):
+            r_start, r_end, c_start, c_end = sei[col]
+            buff = ras.get_row(row, buff)
+            rbuff[c_start:c_end] = buff[c_start:c_end]
+        rast.put_row(rbuff)
+
+
+def patch_map(raster, mapset, mset_str, bbox_list, overwrite=False):
+    """Patch raster using a bounding box list to trim the raster."""
+    # Instantiate the RasterRow input objects
+    rast = RasterRow(raster, mapset)
+    rtype = RasterRow(name=raster, mapset=mset_str % (0, 0))
+    rtype.open('r')
+    rast.open('w', mtype=rtype.mtype, overwrite=overwrite)
+    rtype.close()
+    rasts = []
+    for row, rbbox in enumerate(bbox_list):
+        rrasts = []
+        for col in range(len(rbbox)):
+            rrasts.append(RasterRow(name=raster, mapset=mset_str % (row, col)))
+            rrasts[-1].open('r')
+        rasts.append(rrasts)
+        patch_row(rast, rrasts, rbbox)
+
+    for rrast in rasts:
+        for rast_ in rrast:
+            rast_.close()
+    rast.close()

lib/python/pygrass/modules/grid/split.py

@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  2 19:00:15 2013
+
+@author: pietro
+"""
+from grass.pygrass.gis.region import Region
+from grass.pygrass.vector.basic import Bbox
+
+
+def get_bbox(reg, row, col, width, height, overlap):
+    """Return a Bbox"""
+    north = reg.north - (row * height - overlap) * reg.nsres
+    south = reg.north - ((row + 1) * height + overlap) * reg.nsres
+    east = reg.west + ((col + 1) * width + overlap) * reg.ewres
+    west = reg.west + (col * width - overlap) * reg.ewres
+    return Bbox(north=north if north <= reg.north else reg.north,
+                south=south if south >= reg.south else reg.south,
+                east=east if east <= reg.east else reg.east,
+                west=west if west >= reg.west else reg.west,)
+
+
+def split_region_tiles(region=None, width=100, height=100, overlap=0):
+    """Spit a region into a list of Bbox. ::
+
+        >>> reg = Region()
+        >>> reg.north = 1350
+        >>> reg.south = 0
+        >>> reg.nsres = 1
+        >>> reg.east = 1500
+        >>> reg.west = 0
+        >>> reg.ewres = 1
+        >>> reg.cols
+        1500
+        >>> reg.rows
+        1350
+        >>> split_region_tiles(region=reg, width=1000, height=700, overlap=0)
+        [[Bbox(1350.0, 650.0, 1000.0, 0.0), Bbox(1350.0, 650.0, 1500.0, 1000.0)],
+         [Bbox(650.0, 0.0, 1000.0, 0.0), Bbox(650.0, 0.0, 1500.0, 1000.0)]]
+        >>> split_region_tiles(region=reg, width=1000, height=700, overlap=10)
+        [[Bbox(1350.0, 640.0, 1010.0, 0.0), Bbox(1350.0, 640.0, 1500.0, 990.0)],
+         [Bbox(660.0, 0.0, 1010.0, 0.0), Bbox(660.0, 0.0, 1500.0, 990.0)]]
+    """
+    reg = region if region else Region()
+    ncols = (reg.cols + width - 1) // width
+    nrows = (reg.rows + height - 1) // height
+    box_list = []
+    #print reg
+    for row in xrange(nrows):
+        row_list = []
+        for col in xrange(ncols):
+            #print 'c', c, 'r', r
+            row_list.append(get_bbox(reg, row, col, width, height, overlap))
+        box_list.append(row_list)
+    return box_list