Landscape configuration and composition analysis with GRASS GIS 7

Landscape configuration and composition metrics summarize the structure of the landscape in a determined area and time. These metrics  are closely related to the occurring ecological patterns (e.g. species distribution) and processes (e.g. biogeochemical cycle). As a consequence they are often used as predictors in ecological modelling. Thereby, a researcher who wants to test whether or not the landscape structure has an influence on the studied ecological properties (i.e., species richness, ecosystem functioning, conservation status, energy flow, animal movements, species dispersion, etc), need a tool to derive them from a Land Use or Land Cover (LULC) map.
Norway nike air max

Oakley sunglasses For Sale
Fragstat is by far the most common software used to derive landscape metrics. In this post I am going to describe an Open source and cross-platform alternative with comparable characteristics: the r.li modules in GRASS GIS7.
nike air max running shoes
r.li in GRASS GIS7 has a very user friendly Python graphical interface that allows the calculation of 16 different patch or diversity indexes (here a complete list). Furthermore, if you have little skills in C programming you can add your own personalized index.
yellow oakley sunglasses
To calculate landscape indexes we need a base LULC map, such as ESA Globcover dataset. In this post, I will also use the new Gridded Population of the World V4 (GPW4) to correlate population data to the calculated landscape metrics. Indeed, we are going to test the null hypotheses that landscape structure is not influenced by the human population counts in the study area (California state, USA).
nike air max France
Since importing globcover in GRASS is a bit complicated (due to a problem in the original ESA raster file), I will split this exercise in two parts.

How to import GlobcoverV2.3 and GPWV4 in GRASS

Globcover raster file at 10 seconds (~300 at the equator) spatial resolution can be download here:

or using wget:

mkdir ~/globcover
cd ~/globcover
wget http://due.esrin.esa.int/files/Globcover2009_V2.3_Global_.zip
unzip ~/globcover/Globcover2009_V2.3_Global_.zip
rm ~/globcover/Globcover2009_V2.3_Global_.zip

air jordan wallpaper
After having downloaded and unzipped the file, we open GRASS, create a location with WGS84 as datum and lat-long as projection (EPSG: 4326).

Then, we have to fix an annoying problem in the original globcover raster file due to a bad coordinates extension specification. Indeed, the globcover map extends over 360 degrees of longitude and 180 degrees of latitude. GRASS does not like it! To see the problem you can type:

gdalinfo ~/globcover/GLOBCOVER_L4_200901_200912_V2.3_fixed.tif

Among the information displayed we can read:


Origin = (-180.001388888888897,90.001388888888883)

As a workaround, we use gdalwarp in combination with gdal_traslate to cut the raster in a meaningful coordinates extension and compress it in a rather quick way:

gdalwarp -te -179.9 -89.9 179.9 89.9 -of vrt globcover/GLOBCOVER_L4_200901_200912_V2.3.tif globcover/GLOBCOVER_L4_200901_200912_V2.3_fixed.vrt gdal_translate -co compress=LZW globcover/GLOBCOVER_L4_200901_200912_V2.3_fixed.vrt globcover/GLOBCOVER_L4_200901_200912_V2.3_fixed1.vrt

ladies oakley sunglasses
Now we can import the globcover raster file into GRASS to derive the landscape metrics!

r.in.gdal input=~/globcover/GLOBCOVER_L4_200901_200912_V2.3_fixed1.vrt output=globcover

With the gdal transformation we lost the color table, to reassign it we have to paste the values below in a text file:

11.000000 170:240:240
14.000000 255:255:100
20.000000 220:240:100
30.000000 205:205:102
40.000000 0:100:0
50.000000 0:160:0
60.000000 170:200:0
70.000000 000:60:0
90.000000 40:100:0
100.00000 120:130:0
110.00000 140:160:0
120.00000 190:150:0
130.00000 150:100:0
140.00000 255:180:50
150.00000 255:235:175
160.00000 0:120:90
170.00000 0:150:120
180.00000 0:220:130
190.00000 195:20:0
200.00000 255:245:215
210.00000 0:70:200
220.00000 255:255:255
230.00000 0:0:0

Then save it as ~/globcover/color_table.txt and run in GRASS:

r.colors map=globcover rules=~/globcover/color_table.txt

That’s it! Now we can download GPWv4 here: http://www.ciesin.columbia.edu/data/gpw-v4/

or using wget:

# Download GPW4 using wget
mkdir ~/gpw4
cd ~/gpw4
wget http://www.ciesin.columbia.edu/data/gpw-v4/gpw-v4-preliminary-release-2_population_count_2010.zip
unzip ~/gpw4/gpw-v4-preliminary-release-2_population_count_2010.zip
rm ~/gpw4/gpw-v4-preliminary-release-2_population_count_2010.zip
# Import it
r.in.gdal input=~/gpw4/GL_E_ATOTPOPBT_2010_CNTM.tif out=gpw4
# Set a GRASS nice color table
r.colors map=gpw4 rules=population_dens
globalview
GlobcoverV2.3 in GRASS GIS7
gpw4 in GRASSGIS7 with population_dens color table
GPW4 in GRASSGIS7 with population_dens as color table
Derive landscape metrics with r.li modules from globcover raster map

In this section we focus on r.li.* setup modules and how to calculate landscape metrics using moving windows (other approaches are possible as use a vector file to select areas or draw areas). Before performing landscape metrics analysis, remember to set the region on your study area. The region parameters (extension and resolution) will affect all your further analysis in GRASS!

We will restrict our analysis on California state (USA) to speed up the exercise (a worldwide analysis would have been too long). To download California boundaries we will use again wget:

# Download gadm shape file for USA
mkdir ~/gadm
cd ~/gadm
wget http://biogeo.ucdavis.edu/data/gadm2/shp/USA_adm.zip
unzip ~/gadm/USA_adm.zip
rm ~/gadm/USA_adm.zip
# Import it, removing all the island with a surface less than 100,000 km<sup>2</sup>
v.in.ogr input=~/gadm/USA_adm1.shp out=usa min_area=100000000000 --o snap=1e-06
# Extract California state
v.extract input=usa type=area where='NAME_1="California"' output=ca
# Mask on California (setting a big resolution value to speed up the raster conversion)
g.region vect=ca res=0:00:10 -ap
r.mask vector=ca
# New globcover and population raster for California (Setting the region accordingly to the raster map)
g.region vect=ca res=0:00:10 -ap
r.mapcalc expression="globcover_ca=globcover"
g.region vect=ca res=0:00:30 -ap
r.mapcalc expression="gpw4_ca=int(gpw4)"
# Set the region on globcover resolution since it will be the input map for landscape analysis
g.region vect=ca res=0:00:10 -ap

Before calculating the metrics we must create a setup file that will be used as the input information in all the r.li.* modules. The setup file can be created using r.gui.rlisetup that can be run in the GRASS console or found in “Raster->Landscape patch analysis->Set up analysis and analysis framework”.

g.gui.rlisetup
g.gui.rlisetup 1st window

 

Click on “create” and add all the information about your analysis; we will use moving windows in this exercise then we will fill only the name (openiche_globcover) and “Raster map to use to select areas” box, as shown below.

g.gui.rlisetup_1

 

The next pop-up window is about the kind of analysis we want to perform. In this exercise we will use “Moving Windows” and “Use keyboard to enter sampling area” options. Afterwards we have to define the shape and size of the sampling window, a quite common choice is a rectangular 3×3 pixel window, but it largely depends on the scientific question that we try to answer. In our case a 3×3 rectangular window may be optimal since it mimes the resolution of the population raster (~1 km).

g.gui.rlisetup_4

 

Click on Next and you will see a summary window, where all the properties of your configuration file are shown; click on Next and confirm the creation of the setup file.

g.gui.rlisetup_3

If you see the name of your setup file in the r.gui.rlisetup box you likely did everything correctly, otherwise start again!

Now we have the configuration file ready to be used. We will use it to feed the landscape metric modules which have as name pattern “r.li.<metric>”. Hence, if we want to calculate the patch-density index we look for r.li.patchdensity, and if, instead, we would like to calculate the Renyi index we need to run r.li.renyi, and so on!  A detailed explanation of all the landscape metric modules is available on the always updated GRASSGIS manual pages (i.e, http://grass.osgeo.org/grass70/manuals/r.li.edgedensity.html for edge-density index).

Now let’s run the patch density (number of patches divided surface), edge density (edge length divided surface) and Shannon’s index (the degree of the heterogeneity of the landscape matrix) modules using the configuration file that we just created with g.gui.rlisetup.

r.li.patchdensity in=globcover_ca config=openiche_globcover out=glob_ca_pd
r.li.edgedensity in=globcover_ca config=openiche_globcover out=glob_ca_ed
r.li.shannon in=globcover_ca config=openiche_globcover out=glob_ca_sh

patch density, edge density and Shannon's index in California
Patch density (PD), Edge density(ED) and Shannon’s index (SH) in California

 

The pattern for the three indexes looks uniform in the South East California, with low values, while the variability in the North and Central West part is higher.

In GRASS we can use scatterplot between raster maps using “Create bivariate scatterplot of raster maps” tool in”Analyze map” located in the top of “Map Display” Window (more information here).

scatterplot sh~pd
scatterplot sh~pd
scatterplot ed~pd
scatterplot ed~pd

We see a strong linear relation between ED and PD, while SH shows multiple values near 0 (due to the high dominance of a single Land Cover type) where instead PD has from low to high values (therefore the “rare” Land Cove types in these areas can be grouped in one single patch or sparse in multiple small patches).

We only miss to test if a higher human population also means a significant lower value for the three derived landscape metrics. We can test our null hypotheses importing the raster in R, thanks to spgrass R package, fitting linear models between the raster maps. Anyway, to remove the noise due to the different resolution of the landscape metrics and GPW4 raster maps we have to re-sample them at GPW4 resolution with r.resamp.stats before importing them in R.

# Set the region at GPW4 resolution
g.region vect=ca res=0:00:30 -ap
# Resample using median as resampling method
r.resamp.stats input=glob_ca_pd output=glob_ca_pd_res method=median
r.resamp.stats input=glob_ca_ed output=glob_ca_ed_res method=median
r.resamp.stats input=glob_ca_sh output=glob_ca_sh_res method=median

Now we can go in R, import the raster maps and performing linear regression

# Call R in GRASS
R
# Import raster maps
glob_ca_pd<-readRAST("glob_ca_pd_res")
glob_ca_ed<-readRAST("glob_ca_ed_res")
glob_ca_sh<-readRAST("glob_ca_sh_res")
gpw4_ca<-readRAST("gpw4_ca")
# Linear models
lm_pa<-lm(glob_ca_pd$glob_ca_pd_res~gpw4_ca$gpw4_ca)
lm_ed<-lm(glob_ca_ed$glob_ca_pd_res~gpw4_ca$gpw4_ca)
lm_sh<-lm(glob_ca_sh$glob_ca_sh_res~gpw4_ca$gpw4_ca)

 

Models summary (R2):

GPW~PD R2=0.004 p-value=2e-16

GPW~ED R2=0.003 p-value=2e-16

GPW~SH R2=0.005 p-value=2e-16

To conclude, the null hypotheses of no relationship between landscape metrics and human population is rejected. Indeed, we see a significant (but tiny!) linear negative effect of human population on the landscape metrics. At high human population density (i.e. urban areas) the landscape is rather homogeneous and not diverse. However, the very low R2 values are connected to a more complex situation at very low population density, where the landscape can be more or less homogeneous and diverse, depending on its intrinsic characteristics (e.g. continuous forest, agricultural areas, mixed vegetation, etc). This is clear when looking at the scatterplots below, where at very low population counts we observe a very high variability of landscape metrics:

 

scatterplor_final
Scatterplot Landscape matrics vs population counts

Enjoy GRASS!

One thought on “Landscape configuration and composition analysis with GRASS GIS 7

  1. Thanks for your nice article!

    Concerning the correction of the ESA GlobCover map coordinates, here an alternative solution using “gdal_translate” which does not involve resampling:


    # coords are shifted, fix raster map
    # -a_ullr Assign/override the georeferenced bounds of the output file
    # use larger cache and compress result
    # optionally: write to .vrt instead of .tif
    gdal_translate --config GDAL_CACHEMAX 1200 -a_ullr -180 90 180 -65 \
    -co "COMPRESS=LZW" GLOBCOVER_L4_200901_200912_V2.3.tif GLOBCOVER_L4_200901_200912_V2.3_fixed.tif

    # result:
    gdalinfo GLOBCOVER_L4_200901_200912_V2.3_fixed.tif
    ...
    Origin = (-180.000000000000000,90.000000000000000)
    Pixel Size = (0.002777777777778,-0.002777777777778)
    ...
    Corner Coordinates:
    Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
    Lower Left (-180.0000000, -65.0000000) (180d 0' 0.00"W, 65d 0' 0.00"S)
    Upper Right ( 180.0000000, 90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
    Lower Right ( 180.0000000, -65.0000000) (180d 0' 0.00"E, 65d 0' 0.00"S)
    Center ( 0.0000000, 12.5000000) ( 0d 0' 0.01"E, 12d30' 0.00"N)

Leave a Reply

Your email address will not be published. Required fields are marked *