Analysing the ECA&D climatic data

We'll now explore the prepared dataset and perform some analysis on various maps.

If not yet done, download the provided course data as indicated in the introduction. Copy the given GRASS GIS location to this path: $HOME/grassdata/. The provided location ecad17_ll contains two mapsets:

  1. PERMANENT - mapset with metadata and base cartography (not much actually)
  2. ecad17 - mapset in which all the ECAD datasets are stored
Take a look

Having started GRASS GIS as per introduction with location ecad17_ll and the PERMANENT mapset, you will see the menu of GRASS GIS as well as an empty map canvas. Now, using the menu or icons, load the vector map ne_10m_admin_1_states_provinces. Query the map, zoom, colorize it (double click on the legend entry).

That done, we can close GRASS GIS in order to complete this initial micro-session.

Data organization: First steps with an own mapset

Next, we will now create a new mapset user1 in the location ecad17_ll in which we will perform our own analyses (note: keep original data always separated from own data!). To be able to access in this new mapset user1 the already prepared ECA&D data (stored in mapset ecad17), we will add the ecad17 mapset to the so-called mapset search path. For this, we run the g.mapsets command either from the command line or use the GUI menu.

To support reproducible work, there is a Copy button in the GUI of each module in order to be able to copy the current commands for future replication of the workflow. In case you want to get help about the options and flags of the different commands, use the Help button in the GUI or the --help flag, e.g. g.mapsets --help.

Here the workflow for creating an own mapset etc.:

GRASS GIS logoEnter in a terminal window:

# Launch GRASS GIS, -c creates a new mapset "user1" for our own data, --gui brings up the GUI
grass74 -c $HOME/grassdata/ecad17_ll/user1/ --gui
# List all the raster maps in all mapsets in the current search path
g.list type=raster
### ... how many maps do you see?

# List all accessible mapsets (as per current search path)
g.mapsets -p

# List all available mapsets
g.mapsets -l

# Add the mapset "ecad17" to the search path
g.mapsets mapset=ecad17 operation=add

# List all the mapsets in the now modified search path
g.mapsets -p
# List all accessible raster maps
g.list type=raster

The mapset search path applies to vector maps as well.

TIP: To get out of the "--More--" scrolling, type `q`. Otherwise, scroll up and down as usual.

Country admin borders

Next we import - directly from the online data source naturalearthdata.com - the world admin-level0 data:

# with the special GDAL/OGR drivers vsizip and vsicurl we can directly fetch data from a remote source
# so, we import and check/fix topology on the fly
v.import input="/vsizip/vsicurl/https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip" output=country_boundaries

# add some more metadata
v.support country_boundaries comment="Source: http://www.naturalearthdata.com/downloads/110m-cultural-vectors/"
v.support country_boundaries map_name="Admin0 boundaries from NaturalEarthData.com"

# show attibute table colums
v.info -c country_boundaries
Import of elevation GeoTIFF map

We will now import the elevation map (provided as GeoTIFF files) into the current mapset. Subsequently, we set the color table of elevation map to "elevation" style:

GRASS GIS logoEnter in a terminal window:

# change current directory to where the GeoTIFF files are stored
cd $HOME/geodata/ecad_v17/
# Import the elevation raster map in GeoTIFF format into GRASS GIS
r.in.gdal input=elev_v17.tif output=elev_v17
# Set the color palette of the elevation raster map to 'elevation'
r.colors map=elev_v17 color=elevation
Visualizing maps

Let us now display the ECA&D elevation map and add decorations.

GRASS GIS logoEnter in a terminal window:

# Open a graphical monitor
d.mon wx0

# set current region to elevation map
g.region raster=elev_v17 -p

# tell the GRASS GIS display monitor about the current region change
d.erase

# Display a raster map
d.rast map=elev_v17

# Display a vector map
d.vect map=country_boundaries type=boundary

# Add raster legend (with box around and legend histogram)
d.legend -t -s -b -d raster=elev_v17 title="Elevation [m]" title_fontsize=20 font=sans fontsize=18

# Add north arrow
d.northarrow style=1b text_color=black

# Add text
d.text -b text="Elevation map at ECA&D resolution" color=black bgcolor=229:229:229 align=cc font=sans size=8

TIP: You can use the GUI as well.

ECA&D elevation map
Fig: ECA&D elevation map

Data aggregation in GRASS GIS

In order to dive into some basic data analysis, we next check how to aggregate multiple maps at different time periods based on different statistical methods. Let us create the following aggregated maps from the included monthly long-term averages:

  • annual precipitation for 1951-1980 and 1981-2010
  • mean annual temperature for 1981-2010
  • minimum annual temperature for 1981-2010

GRASS GIS logoEnter in a terminal window:

## precipitation
# Set the computational region to the full raster map (bbox and spatial resolution) and show it (-p print)
g.region raster=precip.1951_1980.01.sum -p
# list all monthly maps for the period 1951-1980 needed for the annual precipitation aggregation
g.list raster pattern="precip.1981_2010.*.sum"
# Now use the r.series command to create annual precip map for the period 1951 to 1980
# note: we now save the output of g.list to a text file and feed it as input to r.series
g.list raster pattern="precip.1981_2010.*.sum" output=list_precip_1981_2010.csv
r.series file=list_precip_1981_2010.csv output=precip.1981_2010.annual.sum method=sum
# Set the color
r.colors precip.1981_2010.annual.sum color=precipitation

# visualize the result
d.mon wx0  # or, d.erase
d.rast precip.1981_2010.annual.sum
d.vect map=country_boundaries type=boundary
d.legend -t -s -b -d raster=precip.1981_2010.annual.sum title="Precipitation [mm]" title_fontsize=20 font=sans fontsize=18
d.text -b text="ECA&D average annual precipitation sums 1981-2010" color=black bgcolor=229:229:229 align=cc font=sans size=8

## mean temperature
# Aggregate the temperature maps average annual temperature
g.list raster pattern="tmean.1981_2010.*.avg"

# note: we now save the output of g.list to a text file and feed it as input to r.series
g.list raster pattern="tmean.1981_2010.*.avg" output=list_tmean_1981_2010.csv
r.series file=list_tmean_1981_2010.csv output=tmean.1981_2010.annual.avg method=average
# Set the color
r.colors tmean.1981_2010.annual.avg color=celsius

# visualize the result
d.erase
d.rast tmean.1981_2010.annual.avg
d.vect map=country_boundaries type=boundary
d.legend -t -s -b -d raster=tmean.1981_2010.annual.avg title="Temperature [deg C]" title_fontsize=20 font=sans fontsize=18
d.text -b text="ECA&D average annual mean temperature 1981-2010" color=black bgcolor=229:229:229 align=cc font=sans size=8

Average annual mean temperatures (1981-2010)
Fig: Average annual mean temperatures (1981-2010)

NOTE: To get rid of the legend or text in the map display canvas, click right-mouse-button on it for a context menu.

Task: Now try to aggregate likewise the mean temperature maps for the period 1951 to 1980 in order to create annual mean temperatures for the same period. Display the annual precipitation and temperature maps with legend and compare the legends using wxGUI or d.mon. For more details on r.series check here.

Univariate statistics in GRASS GIS

Now let's compute univariate statistics of the developed annual maps over the entire area and also over a particular country.

GRASS GIS logoEnter in a terminal window:

# Compute extended univariate statistics
r.univar tmean.1981_2010.annual.avg -e -g
# Compute extended univariate statistics
r.univar precip.1981_2010.annual.sum -e -g
# Get a feel of the vector data and the attributes
v.db.select country_boundaries | head -5
# Find the attributes related to Czech Republic
v.db.select country_boundaries | grep Czech
# Mask out all the area other than Czech Republic
r.mask vector=country_boundaries cats=50

# visualize the difference when having an active MASK
d.erase
d.rast map=elev_v17
d.vect map=country_boundaries type=boundary

# Repeat the univariate statistics, now limited to the Czech Republic
r.univar tmean.1981_2010.annual.avg -e -g
# Repeat the univariate statistics for Czech Republic
r.univar precip.1981_2010.annual.sum -e -g
# Remove the mask
r.mask -r

Univariate statistics of annual precipitation (1981 to 2010)
Fig: Univariate statistics of annual precipitation (1981 to 2010))

Raster map algebra using map calculator

Let us now compute a difference map from the two 30 years aggregated precipitation maps.

GRASS GIS logoEnter in a terminal window:

# we want to be sure to have the computational region set properly       
g.region raster=precip.1981_2010.06.sum -p

# Now use the r.series command to create annual precip map for the period 1951 to 1980
# note: we now save the output of g.list to a text file and feed it as input to r.series
g.list raster pattern="precip.1951_1980.*.sum" output=list_precip_1951_1980.csv
r.series file=list_precip_1951_1980.csv output=precip.1951_1980.annual.sum method=sum
# Set the color
r.colors precip.1951_1980.annual.sum color=precipitation
# Verify stats
r.univar -g precip.1951_1980.annual.sum

# Map algebra using map calculator
r.mapcalc "diff_1951_1981_2010 = precip.1951_1980.annual.sum - precip.1981_2010.annual.sum"
# apply inverse “differences” color table to have dry=red and wet=blue
r.colors map=diff_1951_1981_2010 color=differences -n

# Verify stats, with extended statistics (-e)
r.univar -g -e diff_1951_1981_2010

Differences between precipitation (1951-1980) and (1981-2010)
Fig: Differences between precipitation (1951-1980) and (1981-2010))

Task:Visualize the differences map with the legend and compute univariate statistics limited to the Czech Republic.

Zonal statistics in GRASS GIS

We want to extract a table with basic statistics over all the countries in Europe.

Main steps are:

  • Explore the attributes of the vector data
  • Rasterize the vector data
  • Generate zonal statistics based on r.univar

GRASS GIS logoEnter in a terminal window:

# check the available column (or use wxGUI attribute manager)
v.info -c country_boundaries
# look at country names
v.db.select country_boundaries column="NAME"
# we want to be sure to have the computational region set    
g.region raster=precip.1981_2010.annual.sum -p
# for zonal statistics, we convert the vector polyons to raster model
v.to.rast country_boundaries output=country_boundaries use=cat labelcolumn="NAME"
# check raster polygons
r.category country_boundaries
# compute statistics for each country, store as a CSV table
r.univar -t precip.1981_2010.annual.sum zones=country_boundaries sep="," output=countries_precip_1981_2010_sum.csv
r.univar -t diff_1951_1981_2010 zones=country_boundaries sep="," output=countries_precip_1951_1980_vs_1981_2010_sum.csv

For fun, we can sort the retrieved values (in mm) and generate plot(s).

Column plot of differences between precipitation (1951-1980) and (1981-2010)
Fig: Column plot of differences between precipitation (1951-1980) and (1981-2010))

Up to here we did some basic GIS steps as well as some more advanced processing methods.

Please read on in 03_grass-gis_ecad_randomforest.md.