1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
|
"""
Aggregation methods for space time raster datasets
Usage:
.. code-block:: python
import grass.temporal as tgis
tgis.aggregate_raster_maps(dataset, mapset, inputs, base, start, end, count, method, register_null, dbif)
(C) 2012-2013 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 .datetime_math import create_suffix_from_datetime
from .datetime_math import create_time_suffix
from .datetime_math import create_numeric_suffic
import grass.script as gscript
from grass.exceptions import CalledModuleError
###############################################################################
def collect_map_names(sp, dbif, start, end, sampling):
"""Gather all maps from dataset using a specific sample method
:param sp: The space time raster dataset to select aps from
:param dbif: The temporal database interface to use
:param start: The start time of the sample interval, may be relative or
absolute
:param end: The end time of the sample interval, may be relative or
absolute
:param sampling: The sampling methods to use
"""
use_start = False
use_during = False
use_overlap = False
use_contain = False
use_equal = False
use_follows = False
use_precedes = False
# Initialize the methods
if sampling:
for name in sampling.split(","):
if name == "start":
use_start = True
if name == "during":
use_during = True
if name == "overlap":
use_overlap = True
if name == "contain":
use_contain = True
if name == "equal":
use_equal = True
if name == "follows":
use_follows = True
if name == "precedes":
use_precedes = True
else:
use_start = True
if sp.get_map_time() != "interval":
use_start = True
use_during = False
use_overlap = False
use_contain = False
use_equal = False
use_follows = False
use_precedes = False
where = create_temporal_relation_sql_where_statement(start, end,
use_start,
use_during,
use_overlap,
use_contain,
use_equal,
use_follows,
use_precedes)
rows = sp.get_registered_maps("id", where, "start_time", dbif)
if not rows:
return None
names = []
for row in rows:
names.append(row["id"])
return names
###############################################################################
def aggregate_raster_maps(inputs, base, start, end, count, method,
register_null, dbif, offset=0):
"""Aggregate a list of raster input maps with r.series
:param inputs: The names of the raster maps to be aggregated
:param base: The basename of the new created raster maps
:param start: The start time of the sample interval, may be relative or
absolute
:param end: The end time of the sample interval, may be relative or
absolute
:param count: The number to be attached to the basename of the new
created raster map
:param method: The aggreation method to be used by r.series
:param register_null: If true null maps will be registered in the space
time raster dataset, if false not
:param dbif: The temporal database interface to use
:param offset: Offset to be added to the map counter to create the map ids
"""
msgr = get_tgis_message_interface()
msgr.verbose(_("Aggregating %s raster maps") % (len(inputs)))
output = "%s_%i" % (base, int(offset) + count)
mapset = get_current_mapset()
map_id = output + "@" + mapset
new_map = RasterDataset(map_id)
# Check if new map is in the temporal database
if new_map.is_in_db(dbif):
if gscript.overwrite() is True:
# Remove the existing temporal database entry
new_map.delete(dbif)
new_map = RasterDataset(map_id)
else:
msgr.error(_("Raster map <%(name)s> is already in temporal "
"database, use overwrite flag to overwrite" %
({"name": new_map.get_name()})))
return
msgr.verbose(_("Computing aggregation of maps between %(st)s - %(end)s" % {
'st': str(start), 'end': str(end)}))
# Create the r.series input file
filename = gscript.tempfile(True)
file = open(filename, 'w')
for name in inputs:
string = "%s\n" % (name)
file.write(string)
file.close()
# Run r.series
try:
if len(inputs) > 1000:
gscript.run_command("r.series", flags="z", file=filename,
output=output, overwrite=gscript.overwrite(),
method=method)
else:
gscript.run_command("r.series", file=filename,
output=output, overwrite=gscript.overwrite(),
method=method)
except CalledModuleError:
dbif.close()
msgr.fatal(_("Error occurred in r.series computation"))
# Read the raster map data
new_map.load()
# In case of a null map continue, do not register null maps
if new_map.metadata.get_min() is None and \
new_map.metadata.get_max() is None:
if not register_null:
gscript.run_command("g.remove", flags='f', type='raster',
name=output)
return None
return new_map
##############################################################################
def aggregate_by_topology(granularity_list, granularity, map_list, topo_list,
basename, time_suffix, offset=0, method="average",
nprocs=1, spatial=None, dbif=None, overwrite=False,
file_limit=1000):
"""Aggregate a list of raster input maps with r.series
:param granularity_list: A list of AbstractMapDataset objects.
The temporal extents of the objects are used
to build the spatio-temporal topology with the
map list objects
:param granularity: The granularity of the granularity list
:param map_list: A list of RasterDataset objects that contain the raster
maps that should be aggregated
:param topo_list: A list of strings of topological relations that are
used to select the raster maps for aggregation
:param basename: The basename of the new generated raster maps
:param time_suffix: Use the granularity truncated start time of the
actual granule to create the suffix for the basename
:param offset: Use a numerical offset for suffix generation
(overwritten by time_suffix)
:param method: The aggregation method of r.series (average,min,max, ...)
:param nprocs: The number of processes used for parallel computation
:param spatial: This indicates if the spatial topology is created as
well: spatial can be None (no spatial topology), "2D"
using west, east, south, north or "3D" using west,
east, south, north, bottom, top
:param dbif: The database interface to be used
:param overwrite: Overwrite existing raster maps
:param file_limit: The maximum number of raster map layers that
should be opened at once by r.series
:return: A list of RasterDataset objects that contain the new map names
and the temporal extent for map registration
"""
import grass.pygrass.modules as pymod
import copy
msgr = get_tgis_message_interface()
dbif, connected = init_dbif(dbif)
topo_builder = SpatioTemporalTopologyBuilder()
topo_builder.build(mapsA=granularity_list, mapsB=map_list, spatial=spatial)
# The module queue for parallel execution
process_queue = pymod.ParallelModuleQueue(int(nprocs))
# Dummy process object that will be deep copied
# and be put into the process queue
r_series = pymod.Module("r.series", output="spam", method=[method],
overwrite=overwrite, quiet=True, run_=False,
finish_=False)
g_copy = pymod.Module("g.copy", raster=['spam', 'spamspam'],
quiet=True, run_=False, finish_=False)
output_list = []
count = 0
for granule in granularity_list:
msgr.percent(count, len(granularity_list), 1)
count += 1
aggregation_list = []
if "equal" in topo_list and granule.equal:
for map_layer in granule.equal:
aggregation_list.append(map_layer.get_name())
if "contains" in topo_list and granule.contains:
for map_layer in granule.contains:
aggregation_list.append(map_layer.get_name())
if "during" in topo_list and granule.during:
for map_layer in granule.during:
aggregation_list.append(map_layer.get_name())
if "starts" in topo_list and granule.starts:
for map_layer in granule.starts:
aggregation_list.append(map_layer.get_name())
if "started" in topo_list and granule.started:
for map_layer in granule.started:
aggregation_list.append(map_layer.get_name())
if "finishes" in topo_list and granule.finishes:
for map_layer in granule.finishes:
aggregation_list.append(map_layer.get_name())
if "finished" in topo_list and granule.finished:
for map_layer in granule.finished:
aggregation_list.append(map_layer.get_name())
if "overlaps" in topo_list and granule.overlaps:
for map_layer in granule.overlaps:
aggregation_list.append(map_layer.get_name())
if "overlapped" in topo_list and granule.overlapped:
for map_layer in granule.overlapped:
aggregation_list.append(map_layer.get_name())
if aggregation_list:
msgr.verbose(_("Aggregating %(len)i raster maps from %(start)s to"
" %(end)s") %({"len": len(aggregation_list),
"start": str(granule.temporal_extent.get_start_time()),
"end": str(granule.temporal_extent.get_end_time())}))
if granule.is_time_absolute() is True and time_suffix == 'gran':
suffix = create_suffix_from_datetime(granule.temporal_extent.get_start_time(),
granularity)
output_name = "{ba}_{su}".format(ba=basename, su=suffix)
elif granule.is_time_absolute() is True and time_suffix == 'time':
suffix = create_time_suffix(granule)
output_name = "{ba}_{su}".format(ba=basename, su=suffix)
else:
output_name = create_numeric_suffic(basename, count + int(offset),
time_suffix)
map_layer = RasterDataset("%s@%s" % (output_name,
get_current_mapset()))
map_layer.set_temporal_extent(granule.get_temporal_extent())
if map_layer.map_exists() is True and overwrite is False:
msgr.fatal(_("Unable to perform aggregation. Output raster "
"map <%(name)s> exists and overwrite flag was "
"not set" % ({"name": output_name})))
output_list.append(map_layer)
if len(aggregation_list) > 1:
# Create the r.series input file
filename = gscript.tempfile(True)
file = open(filename, 'w')
for name in aggregation_list:
string = "%s\n" % (name)
file.write(string)
file.close()
mod = copy.deepcopy(r_series)
mod(file=filename, output=output_name)
if len(aggregation_list) > int(file_limit):
msgr.warning(_("The limit of open files (%i) was "\
"reached (%i). The module r.series will "\
"be run with flag z, to avoid open "\
"files limit exceeding."%(int(file_limit),
len(aggregation_list))))
mod(flags="z")
process_queue.put(mod)
else:
mod = copy.deepcopy(g_copy)
mod(raster=[aggregation_list[0], output_name])
process_queue.put(mod)
process_queue.wait()
if connected:
dbif.close()
msgr.percent(1, 1, 1)
return output_list
|