|
@@ -5,7 +5,8 @@
|
|
|
# MODULE: t.rast.neighbors
|
|
|
# AUTHOR(S): Soeren Gebbert
|
|
|
#
|
|
|
-# PURPOSE: Performs a neighborhood analysis for each map in a space time raster dataset.
|
|
|
+# PURPOSE: Performs a neighborhood analysis for each map in a space time
|
|
|
+# raster dataset.
|
|
|
# COPYRIGHT: (C) 2013 by the GRASS Development Team
|
|
|
#
|
|
|
# This program is free software under the GNU General Public
|
|
@@ -71,7 +72,6 @@ import grass.temporal as tgis
|
|
|
|
|
|
############################################################################
|
|
|
|
|
|
-
|
|
|
def main():
|
|
|
|
|
|
# Get the options
|
|
@@ -90,126 +90,101 @@ def main():
|
|
|
dbif = tgis.SQLDatabaseInterfaceConnection()
|
|
|
dbif.connect()
|
|
|
|
|
|
+ overwrite = grass.overwrite()
|
|
|
+
|
|
|
mapset = grass.gisenv()["MAPSET"]
|
|
|
|
|
|
sp = tgis.open_old_space_time_dataset(input, "strds", dbif)
|
|
|
- dummy = sp.get_new_map_instance(None)
|
|
|
-
|
|
|
- temporal_type, semantic_type, title, description = sp.get_initial_values()
|
|
|
- new_sp = tgis.open_new_space_time_dataset(output, "strds", sp.get_temporal_type(),
|
|
|
- title, description, semantic_type,
|
|
|
- dbif, grass.overwrite(), dry=True)
|
|
|
-
|
|
|
- rows = sp.get_registered_maps("id,start_time", where, "start_time", dbif)
|
|
|
+ maps = sp.get_registered_maps_as_objects(where=where, dbif=dbif)
|
|
|
|
|
|
- if not rows:
|
|
|
+ if not maps:
|
|
|
dbif.close()
|
|
|
- grass.fatal(_("Space time raster dataset <%s> is empty") % out_id)
|
|
|
+ grass.fatal(_("Space time raster dataset <%s> is empty") % sp.get_id())
|
|
|
+
|
|
|
+ new_sp = tgis.check_new_space_time_dataset(input, "strds", dbif=dbif,
|
|
|
+ overwrite=overwrite)
|
|
|
|
|
|
count = 0
|
|
|
- proc_count = 0
|
|
|
proc_list = []
|
|
|
|
|
|
- num_rows = len(rows)
|
|
|
- new_maps = {}
|
|
|
+ num_maps = len(maps)
|
|
|
+ new_maps = []
|
|
|
|
|
|
- for row in rows:
|
|
|
+ for map in maps:
|
|
|
count += 1
|
|
|
|
|
|
if count%10 == 0:
|
|
|
- grass.percent(count, num_rows, 1)
|
|
|
+ grass.percent(count, num_maps, 1)
|
|
|
|
|
|
map_name = "%s_%i" % (base, count)
|
|
|
- map_id = dummy.build_id(map_name, mapset)
|
|
|
-
|
|
|
- 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 grass.overwrite():
|
|
|
- # Remove the existing temporal database entry
|
|
|
- new_map.delete(dbif)
|
|
|
- new_map = sp.get_new_map_instance(map_id)
|
|
|
- else:
|
|
|
- grass.error(_("Map <%s> is already in temporal database,"
|
|
|
- " use overwrite flag to overwrite") %
|
|
|
- (new_map.get_map_id()))
|
|
|
- continue
|
|
|
+ new_map = tgis.open_new_map_dataset(map_name, None, mapset,
|
|
|
+ type="raster",
|
|
|
+ temporal_extent=map.get_temporal_extent(),
|
|
|
+ overwrite=overwrite, dbif=dbif)
|
|
|
|
|
|
proc_list.append(Process(target=run_neighbors,
|
|
|
- args=(row["id"],map_name,method,size)))
|
|
|
-
|
|
|
- proc_list[proc_count].start()
|
|
|
- proc_count += 1
|
|
|
+ args=(map.get_id(),map_name,
|
|
|
+ method,size, overwrite)))
|
|
|
+ proc_list[-1].start()
|
|
|
|
|
|
# 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
|
|
|
+ if len(proc_list) == nprocs:
|
|
|
for proc in proc_list:
|
|
|
proc.join()
|
|
|
- exitcodes += proc.exitcode
|
|
|
-
|
|
|
- if exitcodes != 0:
|
|
|
- dbif.close()
|
|
|
- grass.fatal(_("Error while computation"))
|
|
|
+ if proc.exitcode != 0:
|
|
|
+ dbif.close()
|
|
|
+ grass.fatal(_("Error while neighborhood computation"))
|
|
|
|
|
|
# Empty process list
|
|
|
proc_list = []
|
|
|
|
|
|
- # Store the new maps
|
|
|
- new_maps[row["id"]] = new_map
|
|
|
+ # Initlialize and load the content of the map
|
|
|
+ new_maps.append(new_map)
|
|
|
+
|
|
|
+ for proc in proc_list:
|
|
|
+ proc.join()
|
|
|
+ if proc.exitcode != 0:
|
|
|
+ dbif.close()
|
|
|
+ grass.fatal(_("Error while computation"))
|
|
|
|
|
|
- grass.percent(0, num_rows, 1)
|
|
|
+ grass.percent(1, 1, 1)
|
|
|
|
|
|
+ # Open the new space time raster dataset
|
|
|
+ temporal_type, semantic_type, title, description = sp.get_initial_values()
|
|
|
new_sp = tgis.open_new_space_time_dataset(output, "strds",
|
|
|
sp.get_temporal_type(),
|
|
|
title, description,
|
|
|
semantic_type,
|
|
|
- dbif, grass.overwrite(),
|
|
|
- dry=False)
|
|
|
+ dbif, overwrite)
|
|
|
|
|
|
# collect empty maps to remove them
|
|
|
+ num_maps = len(new_maps)
|
|
|
empty_maps = []
|
|
|
|
|
|
# Register the maps in the database
|
|
|
count = 0
|
|
|
- for row in rows:
|
|
|
+ for map in new_maps:
|
|
|
count += 1
|
|
|
+
|
|
|
if count%10 == 0:
|
|
|
- grass.percent(count, num_rows, 1)
|
|
|
- # Register the new maps
|
|
|
- if row["id"] in new_maps:
|
|
|
- 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 new_map.metadata.get_min() is None and \
|
|
|
- new_map.metadata.get_max() is None:
|
|
|
- if not register_null:
|
|
|
- empty_maps.append(new_map)
|
|
|
- continue
|
|
|
-
|
|
|
- old_map = sp.get_new_map_instance(row["id"])
|
|
|
- old_map.select(dbif)
|
|
|
-
|
|
|
- # 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, unit = old_map.get_relative_time()
|
|
|
- new_map.set_relative_time(start, end, unit)
|
|
|
+ grass.percent(count, num_maps, 1)
|
|
|
+
|
|
|
+ # In case of a empty map continue, do not register empty maps
|
|
|
+ map.load()
|
|
|
+ if map.metadata.get_min() is None and \
|
|
|
+ map.metadata.get_max() is None:
|
|
|
+ if not register_null:
|
|
|
+ empty_maps.append(map)
|
|
|
+ continue
|
|
|
|
|
|
- # Insert map in temporal database
|
|
|
- new_map.insert(dbif)
|
|
|
- new_sp.register_map(new_map, dbif)
|
|
|
+ # Insert map in temporal database
|
|
|
+ map.insert(dbif)
|
|
|
+ new_sp.register_map(map, dbif)
|
|
|
|
|
|
# Update the spatio-temporal extent and the metadata table entries
|
|
|
new_sp.update_from_registered_maps(dbif)
|
|
|
- grass.percent(num_rows, num_rows, 1)
|
|
|
+ grass.percent(1, 1, 1)
|
|
|
|
|
|
# Remove empty maps
|
|
|
if len(empty_maps) > 0:
|
|
@@ -224,15 +199,18 @@ def main():
|
|
|
|
|
|
grass.run_command("g.remove", rast=names, quiet=True)
|
|
|
|
|
|
-
|
|
|
dbif.close()
|
|
|
|
|
|
-def run_neighbors(input, output, method, size):
|
|
|
+############################################################################
|
|
|
+
|
|
|
+def run_neighbors(input, output, method, size, overwrite):
|
|
|
"""Helper function to run r.neighbors in parallel"""
|
|
|
return grass.run_command("r.neighbors", input=input, output=output,
|
|
|
- method=method, size=size, overwrite=grass.overwrite(),
|
|
|
+ method=method, size=size,
|
|
|
+ overwrite=overwrite,
|
|
|
quiet=True)
|
|
|
|
|
|
+############################################################################
|
|
|
|
|
|
if __name__ == "__main__":
|
|
|
options, flags = grass.parser()
|