Geografic-Information-System/lib/python/temporal/extract.py

"""!@package grass.temporal

@brief GRASS Python scripting module (temporal GIS functions)

Temporal GIS related functions to be used in Python scripts.

Usage:

@code
import grass.temporal as tgis
@endcode

(C) 2008-2011 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.

@author Soeren Gebbert
"""

from space_time_datasets import *
from multiprocessing import Process

############################################################################

def extract_dataset(input, output, type, where, expression, base, nprocs=1, register_null=False, layer=1, vtype="point,line,boundary,centroid,area,face"):
    """!Extract a subset of a space time raster, raster3d or vector dataset
    
       A mapcalc expression can be provided to process the temporal extracted maps.
       Mapcalc expressions are supported for raster and raster3d maps.
       
       @param input The name of the input space time raster/raster3d dataset 
       @param output The name of the extracted new space time raster/raster3d dataset
       @param type The type of the dataset: "raster", "raster3d" or vector
       @param where The temporal SQL WHERE statement for subset extraction
       @param expression The r(3).mapcalc expression or the v.extract where statement
       @param base The base name of the new created maps in case a mapclac expression is provided 
       @param nprocs The number of parallel processes to be used for mapcalc processing
       @param register_null Set this number True to register empty maps (only raster and raster3d maps)
       @param layer The vector layer number to be used when no timestamped layer is present, default is 1
       @param vtype The feature type to be extracted for vector maps, default is point,line,boundary,centroid,area and face
    """

    # Check the parameters

    if expression and not base:
        core.fatal(_("You need to specify the base name of new created maps"))
    
    mapset =  core.gisenv()["MAPSET"]

    if input.find("@") >= 0:
        id = input
    else:
        id = input + "@" + mapset

    if type == "raster":
	sp = space_time_raster_dataset(id)
    elif type == "raster3d":
	sp = space_time_raster3d_dataset(id)
    elif type == "vector":
	sp = space_time_vector_dataset(id)
	
    dummy = sp.get_new_map_instance(None)
	
    dbif = ()
    dbif.connect()
    
    if sp.is_in_db(dbif) == False:
	dbif.close()
        core.fatal(_("Space time %s dataset <%s> not found") % (type, id))

    if expression and not base:
	dbif.close()
        core.fatal(_("Please specify base="))

    sp.select(dbif)

    if output.find("@") >= 0:
        out_id = output
    else:
        out_id = output + "@" + mapset

    # The new space time dataset
    new_sp = sp.get_new_instance(out_id)
	
    if new_sp.is_in_db():
        if core.overwrite() == False:
	    dbif.close()
            core.fatal(_("Space time %s dataset <%s> is already in database, use overwrite flag to overwrite") % (type, out_id))
    if type == "vector":
	rows = sp.get_registered_maps("id,name,mapset,layer", where, "start_time", dbif)
    else:
	rows = sp.get_registered_maps("id", where, "start_time", dbif)

    new_maps = {}
    if rows:
	num_rows = len(rows)
	
	core.percent(0, num_rows, 1)
	
	# Run the mapcalc expression
        if expression:
	    count = 0
	    proc_count = 0
	    proc_list = []
	    
	    for row in rows:
		count += 1
		
		core.percent(count, num_rows, 1)
		
		map_name = "%s_%i" % (base, count)
		
		# We need to modify the r(3).mapcalc expression
		if type != "vector":
		    expr = "%s = %s" % (map_name, expression)
		    
		    expr = expr.replace(sp.base.get_map_id(), row["id"])
		    expr = expr.replace(sp.base.get_name(), row["id"])
		    
		    # We need to build the id
		    map_id = dummy.build_id(map_name, mapset)
		else:
		    map_id = dummy.build_id(map_name, mapset, row["layer"])

		new_map = sp.get_new_map_instance(map_id)

		# Check if new map is in the temporal database
		if new_map.is_in_db(dbif):
		    if core.overwrite() == True:
			# Remove the existing temporal database entry
			new_map.delete(dbif)
			new_map = sp.get_new_map_instance(map_id)
		    else:
			core.error(_("Map <%s> is already in temporal database, use overwrite flag to overwrite")%(new_map.get_map_id()))
			continue
		
		# Add process to the process list
		if type == "raster":
		    core.verbose(_("Apply r.mapcalc expression: \"%s\"") % expr)
		    proc_list.append(Process(target=run_mapcalc2d, args=(expr,)))
		elif type == "raster3d":
		    core.verbose(_("Apply r3.mapcalc expression: \"%s\"") % expr)
		    proc_list.append(Process(target=run_mapcalc3d, args=(expr,)))
		elif type == "vector":
		    core.verbose(_("Apply v.extract where statement: \"%s\"") % expression)
		    if row["layer"]:
			proc_list.append(Process(target=run_vector_extraction, args=(row["name"] + "@" + row["mapset"], map_name, row["layer"], vtype, expression)))
		    else:
			proc_list.append(Process(target=run_vector_extraction, args=(row["name"] + "@" + row["mapset"], map_name, layer, vtype, expression)))
		
		proc_list[proc_count].start()
		proc_count += 1
		
		# Join processes if the maximum number of processes are reached or the end of the
		# loop is reached
		if proc_count == nprocs or proc_count == num_rows:
		    proc_count = 0
		    exitcodes = 0
		    for proc in proc_list:
			proc.join()
			exitcodes += proc.exitcode
			
		    if exitcodes != 0:
			dbif.close()
			core.fatal(_("Error while computation"))
			
		    # Empty process list
		    proc_list = []
		    
		# Store the new maps
		new_maps[row["id"]] = new_map
	
	core.percent(0, num_rows, 1)
	
	# Insert the new space time dataset
	if new_sp.is_in_db(dbif):
	    if core.overwrite() == True:
		new_sp.delete(dbif)
		new_sp = sp.get_new_instance(out_id)

	temporal_type, semantic_type, title, description = sp.get_initial_values()
	new_sp.set_initial_values(temporal_type, semantic_type, title, description)
	new_sp.insert(dbif)
	
	# collect empty maps to remove them
	empty_maps = []
	
	# Register the maps in the database
        count = 0
        for row in rows:
            count += 1
	    
	    core.percent(count, num_rows, 1)

            old_map = sp.get_new_map_instance(row["id"])
            old_map.select(dbif)
            
            if expression:
		# Register the new maps
		if new_maps.has_key(row["id"]):
		    new_map = new_maps[row["id"]]

		    # Read the raster map data
		    new_map.load()
		    
		    # In case of a empty map continue, do not register empty maps
		    if type == "raster" or type == "raster3d":
			if new_map.metadata.get_min() == None and new_map.metadata.get_max() == None:
			    if not register_null:
				empty_maps.append(new_map)
				continue
		    elif type == "vector":
			if new_map.metadata.get_primitives() == 0 or new_map.metadata.get_primitives() == None:
			    if not register_null:
				empty_maps.append(new_map)
				continue

		    # Set the time stamp
		    if old_map.is_time_absolute():
			start, end, tz = old_map.get_absolute_time()
			new_map.set_absolute_time(start, end, tz)
		    else:
			start, end = old_map.get_relative_time()
			new_map.set_relative_time(start, end)

		    # Insert map in temporal database
		    new_map.insert(dbif)

		    new_sp.register_map(new_map, dbif)
	    else:
		new_sp.register_map(old_map, dbif)          
                
        # Update the spatio-temporal extent and the metadata table entries
        new_sp.update_from_registered_maps(dbif)
	
	core.percent(num_rows, num_rows, 1)
	
	# Remove empty maps
	if len(empty_maps) > 0:
	    names = ""
	    count = 0
	    for map in empty_maps:
		if count == 0:
		    names += "%s"%(map.get_name())
		else:
		    names += ",%s"%(map.get_name())
		count += 1
	    if type == "raster":
		core.run_command("g.remove", rast=names, quiet=True)
	    elif type == "raster3d":
		core.run_command("g.remove", rast3d=names, quiet=True)
	    elif type == "vector":
		core.run_command("g.remove", vect=names, quiet=True)
        
    dbif.close()

###############################################################################

def run_mapcalc2d(expr):
    """Helper function to run r.mapcalc in parallel"""
    return core.run_command("r.mapcalc", expression=expr, overwrite=core.overwrite(), quiet=True)


def run_mapcalc3d(expr):
    """Helper function to run r3.mapcalc in parallel"""
    return core.run_command("r3.mapcalc", expression=expr, overwrite=core.overwrite(), quiet=True)
    

def run_vector_extraction(input, output, layer, type, where):
    """Helper function to run r.mapcalc in parallel"""
    return core.run_command("v.extract", input=input, output=output, layer=layer, type=type, where=where, overwrite=core.overwrite(), quiet=True)