Author: Paulo van Breugel
Updated on: 12/24/2017
About this tutorial
In this tutorial we will look at processing of temporal data, using the aggregated monthly MODIS LST data layers by Metz et al. (2017) as input. We will import the layers in a GRASS GIS database, create a temporal database, and then check out a few options available in GRASS GIS to explore the data. For a more detailed tutorial, see the tutorial Raster time series in GRASS GIS: MODIS Land Surface Temperature (Andrea et al. 2017).
You are expected to be familiar with the basics of data handling and analysis in GRASS GIS. If you are not, please first read Introduction to GRASS GIS and guided tour about working with GRASS GIS. In the tutorials various examples are worked out using the command line (bash) or using a Python script. The tutorial will therefore be easier to follow if you have some knowledge of scripting (bash or Python). If not, see the page on scripting in GRASS GIS on the GRASS GIS Wiki. As written there, GRASS GIS offers two Python Application Programming Interfaces (API), viz. the grass scripting interface and the pygrass api. In this tutorial, I will refer to the former as python and to the latter as pygrass
Metz et al. (2017) published a new dataset based on the land surface temperature (LST) time series based on the MODIS Land Surface Temperature (LST) collection 6. The original data, being remotely sensed data in the thermal range, shows gaps in cloud-covered areas. Metz et al. (2017) present a novel method to fully reconstruct MODIS daily LST products for central Europe at 1 km resolution and globally, at 3 arc-min. The reconstruction was done using the open source tools gdal and GRASS GIS as detailed in the above-mentioned article. The resulting aggregated monthly LST data layers for the period 2003-2016, at a spatial resolution of 3 arc-min, are available for download as geotif files at Zenodo.
Download and import the data
There are different ways to download and import the data. In the example below, I show how to download the files and import them in one go using a small Python script. But, before running the script I had to create the Location with the spatial reference system (SRS) that matches that of the downloaded data.
There are different ways to create the new Location with matching SRS. Use your own preferred way, or check out the example below, in which I used the EPSG code (EPSG 4008) to define the SRS of the new location.
Now we have the location, we can import the data into the location. First step is to download the data from https://zenodo.org/record/1115666/ as zip files (one zipfile per year). You can of course download it manually from the website, but using a script is easier. Try it out yourself using your favorite scripting language, or check out the example Python script below.
Set region and color tables
To check if have imported all layers, you can use the g.list() command or use the convenient Data catalog. Before we continue, note that we did not set the region when creating the Location. So let's set that first to ensure the extent and resolution match that of the imported layers. Or, if you are only interested in a particular area, you can set it to match the bounds of that area.
While we are doing that, we might as well set one color palette for all layers to ensure a matching color legend (see r.colors). Note that in case several input raster maps are provided the range (min, max) of all maps will be used for color table creation. Hence the created color table will span from the smallest minimum to the largest maximum value of all input raster maps and will be applied to all input raster maps.
Creating a map
It is time for a map. You can of course create one using the main GUI. However, below I'll show how to use the display (d.*) functions to generate an png image with the map. This is especially convenient when non-interactive use is needed (e.g., in a script for automated generation of images).
From spatial to space-time raster set
The setup of the spatial framework in GRASS GIS, and the terminology used, may be a bit confusing. You may therefore want to read this introduction to temporal data processing in GRASS GIS on the GRASS GIS wiki. Or, for a short list of the main terms we will encounter in the next sections, see the list below.
Set up the temporal database and temporal raster dataset
Now that we have imported all layers, we can create a space-time raster data set (STRDS) and import the raster layers. First step is to set up the temporal framework for the current mapset (remember, the temporal framework is mapset specific). If you have done this already you can skip this step. To set up the framework, we need to create and set the connection using t.connect.
Next, we can create the STRDS using t.create. To use this function, we need to specify the type of maps, which are rasters in our case. We furthermore need to define the type of time. Time can be absolute (e.g., 1981-12-23 23:30:00) or relative. The latter is represented by an integer (time span) and a time unit (e.g., seconds, day, months). The raster layers we are going to work with represent absolute time periods (Januari 2003 till December 2016).
Register maps into the STRDS
We now have the space-time dataset (STRDS), but it is still empty. If we want to add the LST raster layers to the dataset, we first must add them to the temporal database (this will assign timestamps to the raster layers). Next, we can register them in the space-time dataset (STRDS) we created in the previous section.
There are different ways to register maps in a space-time dataset. The best method will depend on your data and wether you need/have time instances or time intervals (See here for more details). For an overview of the different options, check the t.register manual and the related map registration wiki page.
In the example below, we will use the t.register function to simultaneously register the raster layers in the temporal database and register them within the space time dataset (STRDS) we created in the previous section.
Let's check again the basic info to see how it looks like and then list the raster maps in our LST_Day_monthly STRDS using the code below:
# Check info t.info LST_monthly
- Metz, M., Andreo, V., & Neteler, M. 2017. A New Fully Gap-Free Time Series of Land Surface Temperature from MODIS LST Data. Remote Sensing 9: 1333. URL: http://www.mdpi.com/2072-4292/9/12/1333
- Andrea, Pareeth and van Breugel, 2017. Hands-on to GIS and Remote Sensing with GRASS GIS at ITC - University of Twente on November 3rd, 2017. ITC and OSGeo.nl. http://scientificpubs.com/grass_workshop_itc/index.html
If you have questions
If you have questions or comments about the text, let me know. You can use this contact form. Please make sure to include the page title ("Import and explore reconstructed MODIS LST data") or page name ("import modis lst").