---

Contributed by Agustin Lobo

Contents of this page is available under

Attribution-Share Alike CC-BY-SA 3.

Data Analysis of Exercise Re-vegetation rate after fire. Are there topographic differences? with R

R has hundreds of packages that perform specific tasks and you will use often.
We need rgdal now. For installing it, use the menu or type:
select a repository and proceed. Then type:
to load the package. Note that another package (sp) gets loaded as a dependency. If  by any reason you do not have it installed, you have to install it as well
and repeat require(rgdal).
Once installed, the packages remain installed and you do not have to re-install them in subsequent sessions, but you must use require() to load them for each session
in which you need them. This is done to save memory: you just load what you need.

We start by importing the data, which is the set of shape files that you produced with QGIS:
where dsn is the directory and layer the name of your shape files

ndvidataSDF is a complex spatial object, according to definitions in package sp (this is why rgdal requires sp):
str() reveals the structure of an R object. Note that max.level=2 in str() is important here to abbreviate the otherwise very long output produced by str() for spatial objects.
The "@" are slots that organize the structure of sub-objects within the spatial object. Note that all files of the shape set become one single spatial object in R, so it has to be a complex object.

While we will not be using it, in order to make sure we have really imported what we intended to import, you can plot the spatial object just using plot() (remember R is object-oriented):

Here we will be using just the values which originally were in the *.dbf file, which are the ones in slot @data now, so we select this data for the rest of the exercise:

head(ndvidata) reveals that the table of the vector shape file we are importing has many more columns (variables) than we actually need:
so we make a backup object (ndvidataori) just in case and select the first 10 columns for ndvidata.

We also observe that ELEMENTID is not a unique ID for each case, so we redefine this variable as:

CLOUDS_SUB, the 5th column in this example, is the variable indicating the presence of clouds, which we can check with summary() and table()
Note we select the variable using either its name or its position

In order to discard those observations having clouds, we first make a second object ndvidata2 and then have two alternatives.
The first one is the simplest:


The second alternative recodes "0" as NA and then filters out cases with NA by making use of na.omit():
Note that the result is not the same: 406 cases through the first method, 401 cases through the second.
If you review the output of summary(ndvidata) you observe that there are 6 cases in which ASPECT is NA also (a result of the algorithm used by QGIS to calculate aspect from the DEM layer).
This is obviously dependent upon the actual data you are using, but, in general, this second method is better for this purpose since you discard any observation with at least one NA in any variable:

Now we are going to define the terrain classes according to (http://www.iiasa.ac.at/Research/LUC/GAEZ/landres/terslopes.htm) :

As SLOPE comes from QGIS calculated in degrees, we create a new variable SLOPEPER as
SLOPEPER = 100*tan(SLOPE)

Note that we multiply the angle by pi/180 to convert into radians, which are the units expected by tan(x) (see help(tan))
In a similar way we also create ASPCLAS as a 4-classes vector of orientation:

and now we are ready to create the vector of combined terrain classes TERRCLAS:

A barplot is a good way of visualizing the output of table:

You can actually have both the graphic and the actual values on it:
While boxplot() is a graphic function, it can be used to save (into x, in this case) the x-axis positions of the bars, which we later use to overlay the number of cases for each terrain class.

In both cases, we observe that, with the actual dataset that we imported from the vector shape layer in this case,  the number of cases is sufficient for statistical analysis for all terrain classes.

We plot NDVI values of 2003 and 2005 against the terrain classes:
par() controls many aspects of the plotting style. Here we are using mfrow to split the graphic window into 1 row and 2 colums

The variable of interest is (NDVI2005-NDVI2003F), which we assume is tracking vegetation recover. In order to check the distributions for each terrain class, we take advantage of package lattice, which you must install if you don't have it yet:

As all distributions approach normality, we can run a simple one-way ANOVA:

The ANOVA states that there are significant differences on NDVI change between 2005 and 2003 among terrain classes
(the probability of being wrong if we reject the null hypothesis (=no differences among classes) is 0.0002027)

A boxplot with notches let us visualize which classes are actually different: if the notches of a given pair boxes' notches do not overlap there is strong evidence that their medians differ (see Chambers et al., 1983, Graphical Methods for Data Analysis. Wadsworth & Brooks/Cole. p. 62). The test "Tukey Honest Significant Differences" let us check the associated probabilities:

According to this, moderate slopes facing N vs. flat terrain and, in particular, steeper slopes facing S vs. N are the ones with significant differences in terms of increase of NDVI.

A complementary way of testing the effect of terrain on the recovery of vegetation after fire as recorded by NDVI, is to regress the increase of NDVI against the cosinus of the angle of incidence of solar radiation over the surface (COSI), which we can calculate as (http://www.isprs.org/congresses/beijing2008/proceedings/6b_pdf/40.pdf):

where szen and saz are the solar zenith and azimuth angles, respectively. A figure representing solar angles can be found in http://www.solsticeamateur.com/AzimuthElevation.jpg . For the actual values, dependent upon the day and time, we take the ones provided for the 2003 image.
The idea behind this approach is that assuming a uniform solar radiation on the top of the atmosphere for this region, the differences on the actual incidence of solar radiation on the surface are essentially due to topographic effects. Note that this approach is not considering shading.

COSI is highest for flat terrain and monotonously decreases as the angle of incidence increases:

The relationship of COSI and the previously calculated terrain classes can be explored through the boxplot
Note that N and W facing slopes have lower values of radiation incidence than flat, E and S facing slopes, and that this difference increases with slope.

Finally, we calculate and plot the linear regression:

According to these results, there is a significant effect of COSI on NDVI increase, but this effect accounts for only a 5% of the total varaibility of NDVI increase. This result is consistent with the previous result of the ANOVA by terrain classes, where only moderate slopes facing N vs. flat terrain and steeper slopes facing S vs. N where showing significant differences on NDVI increase.