Forest patches change detection using GRASS and R

As reported here, GRASS 7RC1 has been recently released. In this post I will show how to interface GRASS GIS 7 and R to study the shift in broad-leaved forest land use patches between 1990 and 2006 in Italy. The final result is shown in the figure below. You can see how there has not been so much forest patch movement during the last 20 years in Italy. This is probably due to the quite stable forest coverage in this country (from 28.1 to 30.0 % coverage of land area according to The World Bank data).
green oakley sunglasses

ray ban sun glasses

Shift in forest patches between 1990 and 2006 in Italy.
Shift in forest patches between 1990 and 2006.

Closeout Oakley sunglasses
To create the shift map, first of all we need to download CORINE land cover raster maps for 1990 and 2006 from Copernicus Website. Furthermore, we use wget to download Italy administrative boundary.

nike air max 95 neon

 # Create a new folder in $HOME
 cd $HOME
 mkdir corine
 cd corine
 # Download Italy administrative bounday using wget
 wget -P $HOME/corine/</pre>
 # Move the downloaded zip files in the working directory
 mv $HOME/Downloads/ .
 mv $HOME/Downloads/ .</pre>
 # Unzip the downloaded files
 unzip -d g100_90
 unzip -d g100_06
 unzip -d ITA_adm

 # Reproject ITA boundary vector file in EPSG 3035 using ogr2ogr
 ogr2ogr -t_srs EPSG:3035 ITA_adm/ITA_adm0_3035.shp ITA_adm/ITA_adm0.shp
<a style="color:white" href="" title="nike air jordan v5 retro">nike air jordan v5 retro</a>
 # Create a new GRASS location and mapset where to load the spatial data
 grass70 -c EPSG:3035 $HOME/grassdata/eu_laea/corine

 # Load the CORINE raster files in GRASS in=g100_90/g100_90.tif out=g100_90 in=g100_06/g100_06.tif out=g100_06
<a style="color:white" href="" title="where to buy oakley sunglasses">where to buy oakley sunglasses</a>
 # Load Italian vector file in GRASS in=ITA_adm/ITA_adm0_3035.shp out=ITA_adm

 # Set and check the working region
 g.region raster=g100_90 -p

 # Play with GRASS tool map swiping to have an overview of the change
 g.gui.mapswipe first=g100_90 second=g100_06 mode=swipe

discount ray ban sunglasses for men
GRASS GIS 7 Map swipe tool

After loading all the spatial data in GRASS GIS 7, we have to cut corine only on forests category in Itay.  To complete this step we use several mask and inverse mask operation.

 # Mask on corine 1990
 r.mask g100_90

 # New mask on Forest classes (3.1) corine 1990
 r.mask raster=g100_90 maskcats="23" --o

 # New map for forest classes corine 1990
 r.mapcalc expression='g100_90_c23=g100_90'

 # New mask on Forest classes (3.1) corine 2006
 r.mask raster=g100_06 maskcats=23 --o
<a style="color:white" href="" title="air max tailwind 4">air max tailwind 4</a>
 # New map for forest classes corine 1990
 r.mapcalc expression='g100_06_c23=g100_06_new'

 # 10 km Buffer around Italy, to avoid to exclude forest shift near the boundary
 v.buffer in=ITA_adm out=ITA_adm_10kmbuffer distance=10000

 # Mask on Italy and calculate new CORINE maps.
 r.mask vector=ITA_adm_10kmbuffer
 r.mapcalc expression='g100_90_c23_ITA=g100_90_c23' --o
 r.mapcalc expression='g100_06_c23_ITA=g100_06_c23' --o

 # remove mask
 r.mask -r

 # Set region on new raster
 g.region rast=g100_90_c23_ITA -p

Now the data are ready to analyse pixel distance and angle of shift for each single pixel.

 # Transform corine 1990 raster in point-vector map. input=g100_90_c23_ITA output=g100_90_c23_ITA type=point column=class

 # Transform corine 2006 raster in point-vector map input=g100_06_c23_ITA output=g100_06_c23_ITA type=point column=class

After creating the two point vector maps, we will use v.distance GRASS GIS module to calculate the shift. This operation will take a long time since it will perform a point-by-point shift analysis (~5550000 points). anyway, don’t desperate, GRASS GIS 7 will make your life easier! 😀

 # Add colmuns to the reference map (CORINE 1990) to store shift distance 
 # and angle for each pixel centroids
 v.db.addcolumn map=g100_90_c23_ITA columns="nn_dist double precision",/
 "nn_angle double precision" --o

 # Calculare and store shift distance for each centroid
 v.distance from=g100_90_c23_ITA from_type=point to=g100_06_c23_ITA /
 to_type=point output=forests_shift_italy / upload=dist,to_angle 
 column= nn_dist,nn_angle to_column=class

Even though the distance shift geometries (line for shift, point for no shift) are stored in the vector; they don’t have a category neither a connected table. We the next command we will add categories, table and upload the table with distance and angle of shift.

 # Add the category to each geometry in the vector file.
 v.category input=forests_shift_italy type=line output=forests_shift_italy_new /
 option=add --o

 # Add the table to the vector
 v.db.addtable map=forests_shift_italy_new columns="length DOUBLE", /
 "azimuth DOUBLE","azimuth_rad DOUBLE"

 # Update the table with the azimuth value (angle) for each shift. map=forests_shift_italy_new type=line option=azimuth /
 columns=azimuth units=degrees

 # Update the table with the distance values in meters for each shift. map=forests_shift_italy_new type=line option=length 
 /columns=length units=meters

Finally we have a vector file with data about the distance and angle of each shifted forest patch nicely stored in the table:

cat |length|azimuth|
3152|100 | 270|
3262|100 | 270|
3263|100 | 180|
3375|100 | 270|

It is time to move in R to calculate some descriptive statistics on the shift and to plot the shift map.

 # Let's move in R!
 # Load required packages
 library(sp) # To load and manage spatial data frame
 library(maptools) # To load shape files
 library(ggplot2) # To make figures
 library(raster) # To handle raster data
 library(spgrass) # To interface GRASS GIS 7 with R
 library(jpeg) # To load jpeg in R
 library(plotrix) # To plot circles
 library(colorRamps) # To access cool color ramps

 # To allow plotting of polylines using spplot

 # Run GRASS command in R to select the right files 
 # to be imported (using spgrass function).

 # Import spatial data from the GRASS mapset.

Now we plot our data to see what is the pattern of forest patches shifting in Italy between 1990-2000. We will use topo.color R palette to show shifts with different angles.

 # In order to plot each shift segment (line) with a color correspondent 
 # to its azimuth value, we have to use a loop that iterates on each 
 # geometry stored in the vector file, assigning a color correspondent 
 # to its angle. The loop stores each line geometry in a list together 
 # with the color value. Afterwords we will use this list as layout 
 # in spplot.
 for(i in 1:length(forests_shift_ita_l$cat))
  sub_line = forests_shift_ita_l[i,]
   if(azi==0) {
   lwd=3,col = topo.colors(360)[azi])
 print(paste("I'm in",forests_shift_ita_l$cat[i],sep=" "))

 # Now the vector with line colors and geometry is ready. 
 # We can plot it and save the result in an object
                  scales = list(draw=T,alternating=2),
                  horizontal=FALSE,factor.levels=c("Forests shift 
                  in Italy 1990-2006"), fg="gray",par.strip.text=
                  list(fontfamily='Times')), strip=FALSE)

 # We save the plot in a tiff compressed file with resolution = 300 dpi
                 antialias = "none",pointsize="12")

We also would like to add a cool legend to our shift map. Indeed, while it is intuitive to understand the length of each shift (proportional to the length of each line), the shift direction may be difficult to interpret. Below, I will propose an alternative way to build a legend to address the direction (angle) of each shifted patch of forest. I will add to the legend also the circular median angle shift observed in Italy.

# Circular legend: We use a pie chart with the same color gradient of the shift map.
# We save it as jpeg
 units="mm",antialias = "none",pointsize="12",bg = "white")
 par(mar=c(1, -1, 1, -1))
 pie(rep(1,360), col=topo.colors(360),labels="",clockwise=T,
         border=F,lwd=4,init.angle=90,edges=100,family="Times"), y=0, radius=0.1, nv = 10000, 
         border = "black", col = "white", lty = 1, lwd = 1)
 title(list("Shift angle in degrees",cex=3),line=-1,family="Times")
# Now we import the figure of the circular legend and add it to the spplot
 circular_legend <- readJPEG("circular_legend.jpeg")
 legend <- list("grid.raster", circular_legend,
          x = 4834000, y = 2430000, width = 400000,
          scales = list(draw=T,alternating=2),colorkey=FALSE,
          strip.left=strip.custom(style = 4, 
 horizontal = FALSE,factor.levels=c("Forests shift in Italy 1990-2006"),

# Save it as png

Leave a Reply

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