From e14f69ba4a7aaac45b12d82960beea17b299bcda Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Mon, 9 Sep 2024 14:32:05 -0400 Subject: [PATCH] review time series part 2 --- time_series/time_series_aggregations.qmd | 241 ++++++++---------- ...me_series_management_and_visualization.qmd | 8 +- 2 files changed, 108 insertions(+), 141 deletions(-) diff --git a/time_series/time_series_aggregations.qmd b/time_series/time_series_aggregations.qmd index ba77e1e..c890259 100644 --- a/time_series/time_series_aggregations.qmd +++ b/time_series/time_series_aggregations.qmd @@ -10,13 +10,13 @@ format: code-copy: true code-fold: false categories: [GRASS GIS, Time series, raster, Intermediate, Python] -engine: knitr +engine: jupyter execute: eval: false jupyter: python3 --- -In this second part of the time series tutorials, we will go through different +In this second time series tutorial, we will go through different ways of performing aggregations. Temporal aggregation is a very common task when working with time series as it allows us to summarize information and find patterns of change both in space and time. @@ -24,15 +24,14 @@ find patterns of change both in space and time. ::: {.callout-note title="Setup"} This tutorial can be run locally or in Google Colab. However, make sure you install GRASS GIS 8.4+, download the sample data and set up your project -as explained in the [first part](time_series_management_and_visualization.qmd) -of these time series tutorials. +as explained in the [first](time_series_management_and_visualization.qmd) +time series tutorial. ::: There are two main tools to do time series aggregations in GRASS GIS: [t.rast.aggregate](https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.html) -and [t.rast.series](https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html). +and [t.rast.series](https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html). We'll demonstrate their usage in the upcoming sections. - ```{python} #| echo: false @@ -40,62 +39,58 @@ import os import sys import subprocess -# GRASS GIS database variables -grassbin = "grass-dev" -grassdata = os.path.join(os.path.expanduser('~'), "grassdata") -project = "eu_laea" -mapset = "italy_LST_daily" - +# Ask GRASS GIS where its Python packages are sys.path.append( - subprocess.check_output([grassbin, "--config", "python_path"], text=True).strip() + subprocess.check_output(["grass", "--config", "python_path"], text=True).strip() ) -``` - -```{python} -#| echo: false - # Import the GRASS GIS packages we need import grass.script as gs import grass.jupyter as gj +path_to_project = "eu_laea/italy_LST_daily" + # Start the GRASS GIS Session -session = gj.init(grassdata, project, mapset); +session = gj.init(path_to_project) ``` -## Aggregation with granularity +## Aggregation with granularity -**Granularity** is the greatest common divisor of the temporal extents +**Granularity** is the greatest common divisor of the temporal extents (and possible gaps) of all maps in a space-time dataset. To perform time -series aggregation with granularity, we use +series aggregation with granularity, we use [t.rast.aggregate](https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.html). This tool allows us to aggregate our time series into larger granularities, i.e., -from hourly to daily, from daily to weekly, monthly, etc. It also permits to +from hourly to daily, from daily to weekly, monthly, etc. It also permits to aggregate with *ad hoc* granularities like 3 minutes, 5 days, 3 months, etc. -Supported aggregate methods include average, minimum, maximum, median, etc. See -the [r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) +Supported aggregate methods include average, minimum, maximum, median, etc. See +the [r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) manual page for details. -If you are working with data representing absolute time, *t.rast.aggregate* will -shift the start date for each aggregation process depending on the provided +If you are working with data representing absolute time, *t.rast.aggregate* will +shift the start date for each aggregation process depending on the provided temporal granularity as follows: -- years: will start at the first of January, hence 14-08-2012 00:01:30 will be shifted to 01-01-2012 00:00:00 -- months: will start at the first day of a month, hence 14-08-2012 will be shifted to 01-08-2012 00:00:00 -- weeks: will start at the first day of a week (Monday), hence 14-08-2012 01:30:30 will be shifted to 13-08-2012 01:00:00 -- days: will start at the first hour of a day, hence 14-08-2012 00:01:30 will be shifted to 14-08-2012 00:00:00 -- hours: will start at the first minute of a hour, hence 14-08-2012 01:30:30 will be shifted to 14-08-2012 01:00:00 -- minutes: will start at the first second of a minute, hence 14-08-2012 01:30:30 will be shifted to 14-08-2012 01:30:00 +- **years:** will start at the first of January, hence 14-08-2012 00:01:30 will be shifted to 01-01-2012 00:00:00 +- **months:** will start at the first day of a month, hence 14-08-2012 will be shifted to 01-08-2012 00:00:00 +- **weeks:** will start at the first day of a week (Monday), hence 14-08-2012 01:30:30 will be shifted to 13-08-2012 01:00:00 +- **days:** will start at the first hour of a day, hence 14-08-2012 00:01:30 will be shifted to 14-08-2012 00:00:00 +- **hours:** will start at the first minute of a hour, hence 14-08-2012 01:30:30 will be shifted to 14-08-2012 01:00:00 +- **minutes:** will start at the first second of a minute, hence 14-08-2012 01:30:30 will be shifted to 14-08-2012 01:30:00 -The specification of the temporal relation between the aggregation intervals and -the raster map layers to be aggregated is always formulated from the aggregation -interval viewpoint. By default, it is set to *contains* to aggregate all maps +The specification of the *temporal relation* between the aggregation intervals and +the raster map layers to be aggregated is always formulated from the aggregation +interval viewpoint. By default, it is set to *contains* relation to aggregate all maps that are temporally located in an aggregation interval. +::: {.callout-note title="Temporal relations"} +Something about temporal relations... +::: + ### Monthly and seasonal averages To demonstrate the basic usage of t.rast.aggregate, let's create monthly and seasonal time series starting from the `lst_daily` time series we created -in the [time series management](time_series_management_and_visualization.qmd) +in the [time series management](time_series_management_and_visualization.qmd) tutorial. ```{python} @@ -115,17 +110,17 @@ gs.run_command("t.rast.aggregate", gs.run_command("t.rast.aggregate", input="lst_daily", method="average", - granularity="3 months", - basename="lst_avg", - output="lst_seasonal", + granularity="3 months", + basename="lst_avg", + output="lst_seasonal", suffix="gran", nprocs=4) ``` -If we would like to follow the so called +If we would like to follow the so called [*meteorological seasons*](https://en.wikipedia.org/wiki/Season#Meteorological), i.e., that spring started with March, then we could have used the *where* option -to shift the starting point of the aggregation period as follows: +to shift the starting point of the aggregation period as follows: `where="start_time >= '2013-03-01 00:00:00'"`. Let's compare the granularities using the timeline plot... @@ -139,9 +134,9 @@ Let's compare the granularities using the timeline plot... and use `grass.jupyter` to create a nice animation of the seasonal time series: ```{python} -lstseries = gj.TimeSeriesMap(use_region=True) +lstseries = gj.TimeSeriesMap() lstseries.add_raster_series("lst_seasonal", fill_gaps=False) -lstseries.d_legend(color="black", at=(10,40,2,6)) +lstseries.d_legend(color="black", at=(5, 50, 2, 6)) lstseries.show() ``` @@ -150,10 +145,10 @@ lstseries.show() lstseries.save("lstseries.gif") ``` -For astronomical seasons, i.e., those defined by solstices and equinoxes, we can use -[t.rast.aggregate.ds](https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.ds.html) -as in -[this example](https://grasswiki.osgeo.org/wiki/Temporal_data_processing/seasonal_aggregation), +For astronomical seasons, i.e., those defined by solstices and equinoxes, we can use +[t.rast.aggregate.ds](https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.ds.html) +as in +[this example](https://grasswiki.osgeo.org/wiki/Temporal_data_processing/seasonal_aggregation), or for an approximate result: ```{python} @@ -161,18 +156,16 @@ gs.run_command("t.rast.aggregate", input="lst_daily", method="average", where="start_time >= '2013-03-21 00:00:00'", - granularity="3 months", - basename="lst_avg", - output="lst_seasonal", + granularity="3 months", + basename="lst_avg", + output="lst_seasonal", suffix="gran", nprocs=4) ``` ```{python} # Check info -gs.read_command("t.info", - input="lst_seasonal", - flags="g") +gs.run_command("t.info", input="lst_seasonal") ``` ```{python} @@ -187,14 +180,14 @@ There's an upcoming addon that will perform seasonal aggregations. See: ### Spring warming -We define spring warming as the velocity with which temperature increases from -winter into spring and we can approximate it as the linear regression slope -among LST values of February, March and April. Let's see how to use +We define spring warming as the velocity with which temperature increases from +winter into spring and we can approximate it as the linear regression slope +among LST values of February, March and April. Let's see how to use *t.rast.aggregate* to estimate yearly spring warming values. ```{python} # Define list of months -months=['{0:02d}'.format(m) for m in range(2,5)] +months = [f"{m:02d}" for m in range(2, 5)] print(months) ``` @@ -216,51 +209,31 @@ gs.run_command("t.rast.list", input="annual_spring_warming") ``` ```{python} -spring_warming = gj.TimeSeriesMap(use_region=True) +spring_warming = gj.TimeSeriesMap() spring_warming.add_raster_series("annual_spring_warming", fill_gaps=False) -spring_warming.d_legend(color="black", at=(10,40,2,6)) +spring_warming.d_legend(color="black", at=(5, 50, 2, 6)) spring_warming.show() ``` -:::{.callout-caution title="Question" collapse="true"} -How would you obtain the average spring warming over the whole time series? - -```{python} -# Average spring warming -gs.run_command("t.rast.series", - input="annual_spring_warming", - output="avg_spring_warming", - method="average") - -# Display raster map with interactive class -spw_map = gj.InteractiveMap(width = 500, use_region=True, tiles="CartoDark") -spw_map.add_raster("avg_spring_warming") -spw_map.add_layer_control(position = "bottomright") -spw_map.show() -``` -::: - - ## Full series aggregation -The previous question (as you might have seen) can be solved using -[t.rast.series](https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html). -This tool allows us to aggregate complete time series (or only selected parts) +Tool [t.rast.series](https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html) +allows us to aggregate complete time series (or only selected parts) with a certain method. Available aggregation methods include average, minimum, -maximum, median, etc. See the -[r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) +maximum, median, etc. See the +[r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) manual page for details. ### Maximum and minimum LST We'll demonstrate how to aggregate the whole `lst_daily` time series to obtain the maximum and minimum LST value for the period 2014-2018. Note that the input -is a STRDS object while the output is a single map in which pixel values +is a STRDS object while the output is a single map in which pixel values correspond to the maximum or minimum value of the time series they represent. ```{python} # Get maximum and minimum LST in the STRDS -methods=["maximum","minimum"] +methods = ["maximum", "minimum"] for m in methods: gs.run_command("t.rast.series", @@ -268,7 +241,7 @@ for m in methods: output=f"lst_{m}", method=m, nprocs=4) - + ``` ```{python} @@ -283,59 +256,55 @@ maximum and minimum LST maps. ```{python} # Plot -lst_map=gj.InteractiveMap(width = 500) +lst_map=gj.InteractiveMap() lst_map.add_raster("lst_minimum") lst_map.add_raster("lst_maximum") -lst_map.add_layer_control(position = "bottomright") lst_map.show() ``` ### Long term aggregations -If we want to know how is the temperature on a "typical" January, February and so +If we want to know temperature on a "typical" January, February and so on, we estimate the so called long-term averages or climatologies for each month. -These are usually computed over 20 or 30 years of data, but for the purpose of -demonstrating the usage of *t.rast.series*, we will use our 5-year time series +These are usually computed over 20 or 30 years of data, but for the purpose of +demonstrating the usage of *t.rast.series*, we will use our 5-year time series here. -We will use something we learnt in the listing and selection examples in the +We will use what we learned in the listing and selection examples in the first time series tutorial, i.e., we need to select all maps for each month to perform the aggregation. So, the input will be the whole time series, but we -will use only some maps to estimate the long term monthly average, minimum and -maximum values. Let's see how to do all that with a simple for cycle looping +will use only some maps to estimate the long term monthly average, minimum and +maximum values. Let's see how to do all that with a simple for loop over months and methods: ```{python} # Estimate long term aggr for all months and methods -months=['{0:02d}'.format(m) for m in range(1,13)] -methods=["average","minimum","maximum"] - -for m in months: - for me in methods: - print(f"Aggregating all {m} maps with {me} method") - gs.run_command("t.rast.series", - input="lst_daily", - method=me, - where=f"strftime('%m', start_time)='{m}'", - output="lst_{}_{}".format(me,m)) +months = [f"{m:02d}" for m in range(1, 13)] +methods = ["average", "minimum", "maximum"] + +for month in months: + print(f"Aggregating all {month} maps") + gs.run_command("t.rast.series", + input="lst_daily", + method=methods, + where=f"strftime('%m', start_time)='{month}'", + output=[f"lst_{method}_{month}" for method in methods]) ``` ```{python} # List newly created maps -map_list = gs.list_grouped(type="raster", - pattern="*{average,minimum,maximum}*")['italy_LST_daily'] -print(map_list) +gs.list_grouped(type="raster", pattern="*_{average,minimum,maximum}_*")['italy_LST_daily'] ``` Let's create an animation of the long term average LST. ```{python} # List of average maps -map_list = gs.list_grouped(type="raster", pattern="*average*")['italy_LST_daily'] +map_list = gs.list_grouped(type="raster", pattern="*_average_*")['italy_LST_daily'] # Animation with SeriesMap class -series = gj.SeriesMap(height = 500) +series = gj.SeriesMap() series.add_rasters(map_list) series.d_barscale() series.show() @@ -348,39 +317,38 @@ Why didn't we list with t.rast.list? ### Bioclimatic variables -Perhaps you have heard of [Worldclim](https://www.worldclim.org/) or -[CHELSA](https://chelsa-climate.org/) bioclimatic variables? Well, these are 19 -variables that represent potentially limiting conditions for species. They -derive from the combination of temperature and precipitation long term -aggregations. +Perhaps you have heard of [Worldclim](https://www.worldclim.org/) or +[CHELSA](https://chelsa-climate.org/) bioclimatic variables? Well, these are 19 +variables that represent potentially limiting conditions for species. They +derive from the combination of temperature and precipitation long term +aggregations. Let's use those that we estimated in the previous example to estimate the bioclimatic variables that include temperature. GRASS GIS has a very nice -extension, +extension, [r.bioclim](https://grass.osgeo.org/grass-stable/manuals/addons/r.bioclim.html), to estimate bioclimatic variables. ```{python} # Install extension -gs.run_command("g.extension", - extension="r.bioclim") +gs.run_command("g.extension", extension="r.bioclim") ``` ```{python} # Get lists of maps needed -tmin=gs.list_grouped(type="raster", pattern="lst_minimum_??")["italy_LST_daily"] -tmax=gs.list_grouped(type="raster", pattern="lst_maximum_??")["italy_LST_daily"] -tavg=gs.list_grouped(type="raster", pattern="lst_average_??")["italy_LST_daily"] +tmin = gs.list_grouped(type="raster", pattern="lst_minimum_*")["italy_LST_daily"] +tmax = gs.list_grouped(type="raster", pattern="lst_maximum_*")["italy_LST_daily"] +tavg = gs.list_grouped(type="raster", pattern="lst_average_*")["italy_LST_daily"] -print(tmin,tmax,tavg) +print(tmin, tmax, tavg) ``` ```{python} # Estimate temperature related bioclimatic variables -gs.run_command("r.bioclim", - tmin=tmin, +gs.run_command("r.bioclim", + tmin=tmin, tmax=tmax, - tavg=tavg, - output="worldclim_") + tavg=tavg, + output="worldclim_") ``` ```{python} @@ -388,25 +356,24 @@ gs.run_command("r.bioclim", gs.list_grouped(type="raster", pattern="worldclim*")["italy_LST_daily"] ``` -Let's have a look at some of the maps we just created +Let's have a look at some of the maps we just created: ```{python} # Display raster map with interactive class -bio_map = gj.InteractiveMap(width = 500, use_region=True) +bio_map = gj.InteractiveMap() bio_map.add_raster("worldclim_bio01") bio_map.add_raster("worldclim_bio02") -bio_map.add_layer_control(position = "bottomright") bio_map.show() ``` -Let's assume a certain insect species can only survive where winter mean +Let's assume a certain insect species can only survive where winter mean temperatures are above 5 degrees. How would you determine suitable areas? We'll see that in the upcoming tutorial, stay tuned! ::: {.callout-tip} -If you are curious, maybe you want to have a look at +If you are curious, maybe you want to have a look at [t.rast.mapcalc](https://grass.osgeo.org/grass-stable/manuals/t.rast.mapcalc.html) -and +and [t.rast.algebra](https://grass.osgeo.org/grass-stable/manuals/t.rast.algebra.html). ::: @@ -415,10 +382,10 @@ and - Gebbert, S., Pebesma, E. 2014. _TGRASS: A temporal GIS for field based environmental modeling._ -Environmental Modelling & Software 53, 1-12. +Environmental Modelling & Software 53, 1-12. [DOI](http://dx.doi.org/10.1016/j.envsoft.2013.11.001). - Gebbert, S., Pebesma, E. 2017. _The GRASS GIS temporal framework._ -International Journal of Geographical Information Science 31, 1273-1292. +International Journal of Geographical Information Science 31, 1273-1292. [DOI](http://dx.doi.org/10.1080/13658816.2017.1306862). - [Temporal data processing](https://grasswiki.osgeo.org/wiki/Temporal_data_processing) wiki page. @@ -426,7 +393,7 @@ International Journal of Geographical Information Science 31, 1273-1292. *** :::{.smaller} -The development of this tutorial was funded by the US -[National Science Foundation (NSF)](https://www.nsf.gov/), +The development of this tutorial was funded by the US +[National Science Foundation (NSF)](https://www.nsf.gov/), award [2303651](https://www.nsf.gov/awardsearch/showAward?AWD_ID=2303651). ::: diff --git a/time_series/time_series_management_and_visualization.qmd b/time_series/time_series_management_and_visualization.qmd index b9c1903..ac634d6 100644 --- a/time_series/time_series_management_and_visualization.qmd +++ b/time_series/time_series_management_and_visualization.qmd @@ -10,7 +10,7 @@ format: code-copy: true code-fold: false categories: [GRASS GIS, Time series, raster, Basic, Python] -engine: knitr +engine: jupyter execute: eval: false jupyter: python3 @@ -61,9 +61,9 @@ scheme. In this way, we have: To run this tutorial locally or in Google Colab, you should install GRASS GIS 8.4+, and download the [daily MODIS LST project](https://zenodo.org/records/3564515). -This project contains average daily MODIS LST data reconstructed by -[mundialis GmbH & Co. KG](https://www.mundialis.de/en/) based on -Metz et al (2017). +This project contains average daily MODIS LST (Land Surface Temperature) data +reconstructed by [mundialis GmbH & Co. KG](https://www.mundialis.de/en/) +based on Metz et al (2017). Follow the [Fast track to GRASS GIS]() and [GRASS in Colab]() tutorials to get you started.