-
Notifications
You must be signed in to change notification settings - Fork 170
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
t.rast.climatologies: new addon for GRASS 8 #1005
base: grass8
Are you sure you want to change the base?
Changes from 7 commits
66ac22e
ab97ab7
fc33180
2b855e9
87e03ad
86d57e5
8bfeba5
cfabf7e
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,7 @@ | ||
MODULE_TOPDIR = ../.. | ||
|
||
PGM = t.rast.climatologies | ||
|
||
include $(MODULE_TOPDIR)/include/Make/Script.make | ||
|
||
default: script |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,40 @@ | ||
<h2>DESCRIPTION</h2> | ||
|
||
<b>t.rast.climatologies</b> is a module to calculate climatologies, i.e., long term stats, for single days | ||
or months in a time series. If the <em>s</em> flag is used, the module outputs a new space time | ||
raster dataset of relative temporal type with the aggregated maps. | ||
|
||
<h2>EXAMPLE</h2> | ||
|
||
<h3>Daily climatologies</h3> | ||
|
||
Starting from a space time raster dataset of daily granularity (or granularity lower than one day), | ||
the module will create a new space time raster dataset of relative temporal type containing the | ||
long term average for each day along years. | ||
|
||
<div class="code"> | ||
<pre> | ||
t.rast.climatologies input=myinput output=dailyoutput | ||
</pre> | ||
</div> | ||
|
||
<h3>Monthly climatologies</h3> | ||
|
||
Starting from a space time raster dataset of monthly granularity (or lower than one month), | ||
the module will create two new space time raster datasets containing the long term minimum | ||
and maximum for each month along years. | ||
<div class="code"> | ||
<pre> | ||
t.rast.climatologies input=myinput granularity=month method=min,max output=monthlyoutputmin,monthlyoutputmax | ||
</pre> | ||
</div> | ||
|
||
<h2>SEE ALSO</h2> | ||
|
||
<a href="https://grass.osgeo.org/grass-stable/manuals/r.series.html">r.series</a>, | ||
<a href="https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html">t.rast.series</a> | ||
|
||
|
||
<h2>AUTHOR</h2> | ||
|
||
Luca Delucchi, <i>Fondazione Edmund Mach</i> |
Original file line number | Diff line number | Diff line change | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
@@ -0,0 +1,244 @@ | ||||||||||
#!/usr/bin/env python | ||||||||||
|
||||||||||
################################################ | ||||||||||
# | ||||||||||
# MODULE: t.rast.climatologies | ||||||||||
# AUTHOR(S): Luca Delucchi, Fondazione Edmund Mach | ||||||||||
# PURPOSE: t.rast.climatologies calculates climatologies for space time raster datasets | ||||||||||
# | ||||||||||
# COPYRIGHT: (C) 2023 by Luca Delucchi | ||||||||||
# | ||||||||||
# This program is free software under the GNU General Public | ||||||||||
# License (>=v2). Read the file COPYING that comes with GRASS | ||||||||||
# for details. | ||||||||||
# | ||||||||||
################################################ | ||||||||||
|
||||||||||
# %module | ||||||||||
# % description: Calculates climatologies from a space time raster dataset of absolute temporal type | ||||||||||
# % keyword: temporal | ||||||||||
# % keyword: raster | ||||||||||
# % keyword: aggregation | ||||||||||
# % keyword: series | ||||||||||
# % keyword: time | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option G_OPT_STRDS_INPUT | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option G_OPT_STRDS_OUTPUT | ||||||||||
# % required: no | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option | ||||||||||
# % key: basename | ||||||||||
# % type: string | ||||||||||
# % label: Basename of the newly generated output maps | ||||||||||
# % description: Either a numerical suffix or the start time (s-flag) separated by an underscore will be attached to create a unique identifier | ||||||||||
# % required: yes | ||||||||||
# % multiple: no | ||||||||||
# % gisprompt: | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option | ||||||||||
# % key: method | ||||||||||
# % type: string | ||||||||||
# % description: Aggregate operation to be performed on the raster maps | ||||||||||
# % required: yes | ||||||||||
# % multiple: yes | ||||||||||
# % options: average,count,median,mode,minimum,min_raster,maximum,max_raster,stddev,range,sum,variance,diversity,slope,offset,detcoeff,quart1,quart3,perc90,quantile,skewness,kurtosis | ||||||||||
# % answer: average | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option | ||||||||||
# % key: quantile | ||||||||||
# % type: double | ||||||||||
# % description: Quantile to calculate for method=quantile | ||||||||||
# % required: no | ||||||||||
# % multiple: yes | ||||||||||
# % options: 0.0-1.0 | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option | ||||||||||
# % key: granularity | ||||||||||
# % type: string | ||||||||||
# % label: Aggregate by day or month | ||||||||||
# % required: yes | ||||||||||
# % multiple: no | ||||||||||
# % options: day, month | ||||||||||
# % answer: day | ||||||||||
# %end | ||||||||||
|
||||||||||
# %option | ||||||||||
# % key: nprocs | ||||||||||
# % type: integer | ||||||||||
# % description: Number of processes to run in parallel | ||||||||||
# % required: no | ||||||||||
# % multiple: no | ||||||||||
# % answer: 1 | ||||||||||
# %end | ||||||||||
|
||||||||||
# %flag | ||||||||||
# % key: s | ||||||||||
# % description: Do not create a space time raster dataset as output | ||||||||||
# %end | ||||||||||
|
||||||||||
import copy | ||||||||||
from datetime import datetime | ||||||||||
|
||||||||||
|
||||||||||
def main(): | ||||||||||
import grass.pygrass.modules as pymod | ||||||||||
import grass.script as gscript | ||||||||||
import grass.temporal as tgis | ||||||||||
|
||||||||||
options, flags = gscript.parser() | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
strds = options["input"] | ||||||||||
output = options["output"] | ||||||||||
method = options["method"] | ||||||||||
gran = options["granularity"] | ||||||||||
basename = options["basename"] | ||||||||||
nprocs = options["nprocs"] | ||||||||||
quantile = options["quantile"] | ||||||||||
space_time = flags["s"] | ||||||||||
|
||||||||||
# check if quantile value is used correctly | ||||||||||
if "quantile" in method and not quantile: | ||||||||||
gscript.fatal(_("Number requested methods and output maps do not match.")) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
elif quantile and "quantile" not in method: | ||||||||||
gscript.warning( | ||||||||||
_("Quantile option set but quantile not selected in method option") | ||||||||||
) | ||||||||||
Comment on lines
+109
to
+111
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
|
||||||||||
# Check if number of methods and output maps matches | ||||||||||
if "quantile" in method: | ||||||||||
len_method = len(method.split(",")) - 1 | ||||||||||
else: | ||||||||||
len_method = len(method.split(",")) | ||||||||||
|
||||||||||
if not space_time: | ||||||||||
if (len(list(filter(None, quantile.split(",")))) + len_method) != len( | ||||||||||
output.split(",") | ||||||||||
): | ||||||||||
gscript.fatal(_("Number requested methods and output maps do not match.")) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
|
||||||||||
tgis.init() | ||||||||||
# We need a database interface | ||||||||||
dbif = tgis.SQLDatabaseInterfaceConnection() | ||||||||||
dbif.connect() | ||||||||||
|
||||||||||
insp = tgis.open_old_stds(strds, "strds", dbif) | ||||||||||
temporal_type, semantic_type, title, description = insp.get_initial_values() | ||||||||||
if temporal_type != "absolute": | ||||||||||
gscript.fatal( | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
_( | ||||||||||
"Space time raster dataset temporal type is not absolute, this module requires absolute time" | ||||||||||
) | ||||||||||
) | ||||||||||
maps = insp.get_registered_maps_as_objects(None, "start_time", dbif) | ||||||||||
if maps is None: | ||||||||||
gscript.fatal( | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
_( | ||||||||||
"No maps selected in the space time raster dataset {};" | ||||||||||
" it might be empty or the where option returns no data".format(strds) | ||||||||||
) | ||||||||||
) | ||||||||||
# start the r.series module to be used in a ParallelModuleQueue | ||||||||||
mod = pymod.Module("r.series") | ||||||||||
mod.inputs.method = method | ||||||||||
mod.flags.quiet = True | ||||||||||
if quantile: | ||||||||||
mod.inputs.quantile = quantile | ||||||||||
process_queue = pymod.ParallelModuleQueue(int(nprocs)) | ||||||||||
mapset = tgis.core.get_current_mapset() | ||||||||||
# depending on granularity it calculates daily or monthly climatologies | ||||||||||
outmaps = [] | ||||||||||
if gran == "day": | ||||||||||
outunit = "days" | ||||||||||
# for each day | ||||||||||
for doy in range(1, 367): | ||||||||||
doystr = datetime.strptime(f"2000 {doy}", "%Y %j").strftime("%m_%d") | ||||||||||
thiswhere = f"strftime('%m_%d', start_time) == '{doystr}'" | ||||||||||
selemaps = insp.get_registered_maps_as_objects( | ||||||||||
thiswhere, "start_time", dbif | ||||||||||
) | ||||||||||
maps_name = [sam.get_id() for sam in selemaps] | ||||||||||
# check if there are maps for that day | ||||||||||
if len(maps_name) > 0: | ||||||||||
outname = f"{basename}_{doystr}" | ||||||||||
runmod = copy.deepcopy(mod) | ||||||||||
runmod.inputs.input = ",".join(maps_name) | ||||||||||
runmod.outputs.output = outname | ||||||||||
process_queue.put(runmod) | ||||||||||
map_layer = tgis.space_time_datasets.RasterDataset( | ||||||||||
f"{outname}@{mapset}" | ||||||||||
) | ||||||||||
extent = tgis.RelativeTemporalExtent( | ||||||||||
start_time=doy - 1, | ||||||||||
end_time=doy, | ||||||||||
unit=outunit, | ||||||||||
) | ||||||||||
map_layer.set_temporal_extent(extent=extent) | ||||||||||
outmaps.append(map_layer) | ||||||||||
|
||||||||||
if doy % 10 == 0: | ||||||||||
gscript.percent(doy, 366, 1) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
else: | ||||||||||
outunit = "months" | ||||||||||
for month in range(1, 13): | ||||||||||
monthstr = "{:02d}".format(month) | ||||||||||
thiswhere = f"strftime('%m', start_time) == '{monthstr}'" | ||||||||||
selemaps = insp.get_registered_maps_as_objects( | ||||||||||
thiswhere, "start_time", None | ||||||||||
) | ||||||||||
maps_name = [sam.get_id() for sam in selemaps] | ||||||||||
if len(maps_name) > 0: | ||||||||||
outname = f"{basename}_{monthstr}" | ||||||||||
runmod = copy.deepcopy(mod) | ||||||||||
runmod.inputs.input = ",".join(maps_name) | ||||||||||
runmod.outputs.output = outname | ||||||||||
process_queue.put(runmod) | ||||||||||
map_layer = tgis.space_time_datasets.RasterDataset( | ||||||||||
f"{outname}@{mapset}" | ||||||||||
) | ||||||||||
extent = tgis.RelativeTemporalExtent( | ||||||||||
start_time=month - 1, | ||||||||||
end_time=month, | ||||||||||
unit=outunit, | ||||||||||
) | ||||||||||
map_layer.set_temporal_extent(extent=extent) | ||||||||||
outmaps.append(map_layer) | ||||||||||
gscript.percent(month, 12, 1) | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
|
||||||||||
if not space_time: | ||||||||||
# create new space time raster dataset | ||||||||||
output_strds = tgis.open_new_stds( | ||||||||||
output, | ||||||||||
"strds", | ||||||||||
"relative", | ||||||||||
f"{gran} {method} climatologies", | ||||||||||
f"Climatologies created from {strds}, {gran} {method} maps", | ||||||||||
semantic_type, | ||||||||||
dbif, | ||||||||||
gscript.overwrite(), | ||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. [ruff] reported by reviewdog 🐶
Suggested change
|
||||||||||
) | ||||||||||
register_null = False | ||||||||||
# register maps into space time raster dataset | ||||||||||
tgis.register_map_object_list( | ||||||||||
"rast", | ||||||||||
outmaps, | ||||||||||
output_strds, | ||||||||||
register_null, | ||||||||||
outunit, | ||||||||||
dbif, | ||||||||||
) | ||||||||||
|
||||||||||
# Update the raster metadata table entries | ||||||||||
output_strds.metadata.update(dbif) | ||||||||||
|
||||||||||
dbif.close() | ||||||||||
return True | ||||||||||
|
||||||||||
|
||||||||||
if __name__ == "__main__": | ||||||||||
main() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
[ruff] reported by reviewdog 🐶