Import and explore reconstructed MODIS LST data

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

The data

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.

You can also import it in a location with different SRS using r.import, which offers on the fly reprojection.

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.

Create a new location with a SRS that matches that of the downloaded layers.

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.

Download the data layers and import them in your GRASS GIS database

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.

First, set the region to match that of the imported layers. Next, set one color palette for all layers

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).

Create a map of the LST of January 2003. Add a legend with histogram and a title at the bottom.


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.

Click to expand / collapse list...

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).

After we have created the STRDS, we can use the t.list and t.info functions to get more information about the space time dataset.

Set up a new temporal framework and create the STDS

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.

Register the LTS maps in the STRDS

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


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").