Skip to content

Commit

Permalink
t.rast.univar / t.rast3d.univar: Add support for zones (#2588)
Browse files Browse the repository at this point in the history 
* add support for zones

* add support for zones

* add test for zones

* clean properly

* add test for zones

* add support for zones

* add credits

* add credits

* add zones in manual

* add check for zones map

* black

* use RasterRow context manager

* fix zones existence test

* try if centos fails because of context manager

* fix indent

* change import order

* zones check for t.rast3d.univar

* avoid RasterRow

* avoid RasterRow

* fix if

* fix if

* fix raster_info

* fix raster_info

* fix raster_info

* fix raster_info

* import as gs

* import as gs and zones

* string formating

* clean test

* add semantic labels

* black

* remove gscript

* name tests

* name tests

* move parser, use kwargs

* move parser, use kwargs, allow DCELL

* Import tgis after parser

* Import tgis after parser

* Fix typo
  • Loading branch information
@ninsbl
ninsbl committed Oct 24, 2022
1 parent b2435bf commit 0e6fd7c
Show file tree
Hide file tree
Showing 7 changed files with 349 additions and 78 deletions.
122 changes: 89 additions & 33 deletions python/grass/temporal/univar_statistics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,15 @@


def print_gridded_dataset_univar_statistics(
type, input, output, where, extended, no_header=False, fs="|", rast_region=False
type,
input,
output,
where,
extended,
no_header=False,
fs="|",
rast_region=False,
zones=None,
):
"""Print univariate statistics for a space time raster or raster3d dataset
Expand All @@ -43,8 +51,14 @@ def print_gridded_dataset_univar_statistics(
:param rast_region: If set True ignore the current region settings
and use the raster map regions for univar statistical calculation.
Only available for strds.
:param zones: raster map with zones to calculate statistics for
"""

stats_module = {
"strds": "r.univar",
"str3ds": "r3.univar",
}[type]

# We need a database interface
dbif = SQLDatabaseInterfaceConnection()
dbif.connect()
Expand All @@ -54,9 +68,14 @@ def print_gridded_dataset_univar_statistics(
if output is not None:
out_file = open(output, "w")

rows = sp.get_registered_maps("id,start_time,end_time", where, "start_time", dbif)
strds_cols = (
"id,start_time,end_time,semantic_label"
if type == "strds"
else "id,start_time,end_time"
)
rows = sp.get_registered_maps(strds_cols, where, "start_time", dbif)

if not rows:
if not rows and rows != [""]:
dbif.close()
err = "Space time %(sp)s dataset <%(i)s> is empty"
if where:
Expand All @@ -66,40 +85,60 @@ def print_gridded_dataset_univar_statistics(
)

if no_header is False:
string = ""
string += "id" + fs + "start" + fs + "end" + fs + "mean" + fs
string += "min" + fs + "max" + fs
string += "mean_of_abs" + fs + "stddev" + fs + "variance" + fs
string += "coeff_var" + fs + "sum" + fs + "null_cells" + fs + "cells"
string += fs + "non_null_cells"
cols = (
["id", "semantic_label", "start", "end"]
if type == "strds"
else ["id", "start", "end"]
)
if zones:
cols.append("zone")
cols.extend(
[
"mean",
"min",
"max",
"mean_of_abs",
"stddev",
"variance",
"coeff_var",
"sum",
"null_cells",
"cells",
"non_null_cells",
]
)
if extended is True:
string += fs + "first_quartile" + fs + "median" + fs
string += "third_quartile" + fs + "percentile_90"
cols.extend(["first_quartile", "median", "third_quartile", "percentile_90"])
string = fs.join(cols)

if output is None:
print(string)
else:
out_file.write(string + "\n")

flag = "g"

if extended is True:
flag += "e"
if type == "strds" and rast_region is True:
flag += "r"

for row in rows:
string = ""
id = row["id"]
start = row["start_time"]
end = row["end_time"]
semantic_label = (
""
if type != "strds" or not row["semantic_label"]
else row["semantic_label"]
)

flag = "g"

if extended is True:
flag += "e"
if type == "strds" and rast_region is True:
flag += "r"

if type == "strds":
stats = gscript.parse_command("r.univar", map=id, flags=flag)
elif type == "str3ds":
stats = gscript.parse_command("r3.univar", map=id, flags=flag)
univar_stats = gscript.read_command(
stats_module, map=id, flags=flag, zones=zones
).rstrip()

if not stats:
if not univar_stats:
if type == "strds":
gscript.warning(
_("Unable to get statistics for raster map " "<%s>") % id
Expand All @@ -109,19 +148,36 @@ def print_gridded_dataset_univar_statistics(
_("Unable to get statistics for 3d raster map" " <%s>") % id
)
continue
eol = ""

string += str(id) + fs + str(start) + fs + str(end)
string += fs + str(stats["mean"]) + fs + str(stats["min"])
string += fs + str(stats["max"]) + fs + str(stats["mean_of_abs"])
string += fs + str(stats["stddev"]) + fs + str(stats["variance"])
string += fs + str(stats["coeff_var"]) + fs + str(stats["sum"])
string += fs + str(stats["null_cells"]) + fs + str(stats["cells"])
string += fs + str(int(stats["cells"]) - int(stats["null_cells"]))
if extended is True:
string += fs + str(stats["first_quartile"]) + fs + str(stats["median"])
for idx, stats_kv in enumerate(univar_stats.split(";")):
stats = gscript.utils.parse_key_val(stats_kv)
string += (
fs + str(stats["third_quartile"]) + fs + str(stats["percentile_90"])
f"{id}{fs}{semantic_label}{fs}{start}{fs}{end}"
if type == "strds"
else f"{id}{fs}{start}{fs}{end}"
)
if zones:
if idx == 0:
zone = str(stats["zone"])
string = ""
continue
string += f"{fs}{zone}"
if "zone" in stats:
zone = str(stats["zone"])
eol = "\n"
else:
eol = ""
string += f'{fs}{stats["mean"]}{fs}{stats["min"]}'
string += f'{fs}{stats["max"]}{fs}{stats["mean_of_abs"]}'
string += f'{fs}{stats["stddev"]}{fs}{stats["variance"]}'
string += f'{fs}{stats["coeff_var"]}{fs}{stats["sum"]}'
string += f'{fs}{stats["null_cells"]}{fs}{stats["n"]}'
string += f'{fs}{stats["n"]}'
if extended is True:
string += f'{fs}{stats["first_quartile"]}{fs}{stats["median"]}'
string += f'{fs}{stats["third_quartile"]}{fs}{stats["percentile_90"]}'
string += eol

if output is None:
print(string)
Expand Down
9 changes: 8 additions & 1 deletion temporal/t.rast.univar/t.rast.univar.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,15 @@ <h2>DESCRIPTION</h2>
By default it returns the name of the map, the start and end date of
dataset and the following values: mean, minimum and maximum vale,
mean_of_abs, standard deviation, variance, coeff_var, number of null
cells, total number of cell.
cells, total number of cells.
<p>
Using the <em>e</em> flag it can calculate also extended statistics:
first quartile, median value, third quartile and percentile 90.
<p>
If a <em>zones</em> raster map is provided, statistics are computed for
each zone (category) in that input raster map. The <em>zones</em> option
does not support Spatio-Temporal-Raster-Datasets (STRDS) but only a single,
static raster map.

<h2>EXAMPLE</h2>

Expand All @@ -32,8 +37,10 @@ <h2>SEE ALSO</h2>
<em>
<a href="t.create.html">t.create</a>,
<a href="t.info.html">t.info</a>
<a href="r.univar.html">t.create</a>,
</em>

<h2>AUTHOR</h2>

S&ouml;ren Gebbert, Th&uuml;nen Institute of Climate-Smart Agriculture
Stefan Blumentrath, (Support for zones)
33 changes: 27 additions & 6 deletions temporal/t.rast.univar/t.rast.univar.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,12 @@
# %option G_OPT_STRDS_INPUT
# %end

# %option G_OPT_R_INPUT
# % key: zones
# % description: Raster map used for zoning, must be of type CELL
# % required: no
# %end

# %option G_OPT_F_OUTPUT
# % required: no
# %end
Expand Down Expand Up @@ -60,24 +66,27 @@
# % guisection: Formatting
# %end

import grass.script as grass

import grass.script as gs

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


def main():
# Get the options and flags
options, flags = gs.parser()

# lazy imports
import grass.temporal as tgis

# Get the options
# Define variables
input = options["input"]
zones = options["zones"]
output = options["output"]
where = options["where"]
extended = flags["e"]
no_header = flags["u"]
rast_region = bool(flags["r"])
separator = grass.separator(options["separator"])
separator = gs.separator(options["separator"])

# Make sure the temporal database exists
tgis.init()
Expand All @@ -87,11 +96,23 @@ def main():
if output == "-":
output = None

# Check if zones map exists and is of type CELL
if zones:
if gs.raster.raster_info(zones)["datatype"] != "CELL":
gs.fatal(_("Zoning raster must be of type CELL"))

tgis.print_gridded_dataset_univar_statistics(
"strds", input, output, where, extended, no_header, separator, rast_region
"strds",
input,
output,
where,
extended,
no_header=no_header,
fs=separator,
rast_region=rast_region,
zones=zones,
)


if __name__ == "__main__":
options, flags = grass.parser()
main()

0 comments on commit 0e6fd7c

Please sign in to comment.