Automated Mapping in R

After creating our base-layers within a PostGIS database as described here we used R to map our layers and create a plot showing the area count for each year from 1920 to 2012. In total we created 1840 plots (92 years x 20 buffers) that where then combined to produce this video.

The source code for the mappings is complex and wont be discussed in detail here but a brief description of the single steps as well as the code is provided for anybody interested in R-Programming:

A. Session Setup
Here we set our working directory, load all required libraries and suppress scientific notations (later-on important for the plots)

B. Load the data
We load all relevant data.

  • ucs_cum is a table that contains the area counts for all years and all protected area classes. We will need it for the area count on the left side of the video.
  • brasil and world.countries are two shapefiles that we need for our plots. World.countries contains the administration boarders for  all countries worldwide and brasil the boarders for Brazil. We fortify these shapefiles to make them understandable and plottable for ggplot2.
  • map.picture3 is the picture on the upper left of the video. We rasterize it to make it plottable.
  • Next we load all shapefiles that we need for the maps. All in all we have 3 different protected areas with each 20 shapefiles. To avoid too much code we have to automatize this follwing a simple three step approach: First we create a list with all filenames in a loop. Second we create an empty list where we will store our imports. Third we import the shapefiles with a loop thereby using the list of file-names. shapefiles.list.import contains all layers for all multiple use areas. shapefile.list.upi contains all layers for strictly protected areas and shapefile.list.ind contains all layers for indigenous territories. Afterwards we rearrange the order of the list and plot the minimum years because we will need this information later.

C: Construct parameters for the lineplot

Here we calculate the extension of all protected areas (total_protected_area_sqkm) and the extension of indigenous areas + strictly protected areas (total_uus_ind_sqkm). We need this information later  for the areacount plot.

D: Define viewports for the plots

We define different Viewports for the plots. Inside these Viewports we will later paste the different parts of our video (the map, the area count and the picture)

E: Define color schemes for the plots

We define custom colors for the different protected areas.

F: Define the different themes for the plots

We define how our plots made with ggplot2 should be lay-outed. Actually this piece of code was only created after we created the plots which come later. We externalize this part of the code because the plot function that follows afterwards is very long. If you do testplots and want to adjust small parameters in the layout, it might be an advantage to have the layout outside of the plot function.

G: Define additional parameters for the plots

Here we define the plot-title.

H: Construct the plotfunction

Now it gets tricky. In this function we do the following things.

  1. First we create additional layers that contain the information about protected areas (tmp.data.unbuffered.1-3). These layers will display all protected areas that have been created until the respective year. We have to do it with an if-else statement since we don’t know if there is any data for the given year.
  2. Second we create additional layers that contain information only about the respective year. In this layer we plot our buffers that indicate the newly created areas (tmp.data1-3). Again we have to do it with an if-else statement since we don’t know if there is any data for the given year.
  3. We construct the mapplot by adding the baselayers (countries worldwide, Brazil) and additional layers from 1 and 2. We furthermore add limits, title, projection and theme configurations to the plot.
  4. We construct the lineplot (areacount) with geom_ribbon from ggplot2.
  5. We define the area count as a title of the lineplot and modify the plot.
  6. We load the picture
  7. We plot our plots. Therefore we first open a new png file. Inside this file we create a Grid where we push the Viewports. In each Viewport we place our plots.
  8. We turn the device off and finish our plotting function.

I: Plot the data

We plot the data iteratively over all years and all buffersizes. We therefore define the output path and create two vectors for the iteration. One containing the years to be plotted and another one containing the alpha levels for the different buffer sizes (if you watch the video carefully you might notice that the alpha level changes with buffersize).

We then plot everything which takes about 4 hours to complete.

J: Close the session

We close the session and save all created objects from the Global environment.

 

# protected_areas_plotting_1
# Author: Johannes Schielein
# Outline:
  # A. Session Setup
  # B. Load the data
  # C: Construct parameters for the lineplot
  # D: Define viewports for the plots
  # E: Define color schemes for the plots
  # F: Define the different themes for the plots
  # G: Define additional parameters for the plots
  # H: Construct the plotfunctions
  # I: Plot the data
  # J: Close the session

# ----------- A. Session Setup -----------
# set wd
setwd("~/nas/documents/programming/R/protected_areas_plotting/")

# load workspace (if already saved once)
load("r-workspaces/protected_areas_plotting_1.Rdata")

# activate required libraries
library("rgdal")
library("RODBC")
library("mapproj")
library("ggplot2")
library("gridExtra")
library("scales")
library("knitr")
library("extrafont")
#font_import() # only if not already imported
library("png")

# surpress cientific notation
options(scipen=999)

# ----------- B.Load the data -----------
### Load the base data
# load cumulative uc data
ucs_cum<-read.csv("statistics/statistics3112.csv")

# shp brasil
brasil<-readOGR("World Countries 2/","subset_brazil")
brasil.fortified<-fortify(brasil)

# shp worldcountries
world.countries<-readOGR("World Countries 2/","sa")
world.countries.fortified<-fortify(world.countries)

### load the PNG foto for the plots
map.picture3<-readPNG("pics/fotos3.png") #load fotos
map.picture3.rastered <-rasterGrob(map.picture3)

### Load all uus shapefiles
## create an interation vector with the names endings
# the endings of all shapfile names
shapefile.iteration.vector<-vector(length=20) #create the vector
shapefile.iteration.vector[1]<-0 # set the first place zero
for (i in 2:20)shapefile.iteration.vector[i]<-shapefile.iteration.vector[i-1]+5 #set the rest places according to the names endings
shapefile.iteration.vector #show the vector
## creat a vector with the complete names
shapefile.list<-vector(length=20) # create the vector
# create a list with the shapefile names
for (i in shapefile.iteration.vector)shapefile.list[i/5]<-paste("_uus_buffer",i,sep="") # paste beginning and end of names
shapefile.list[20]<-("_uus_buffer0") # _uus_buffer0 was lost so set it automatically
shapefile.list
## import all shapefiles to the list
shapefiles.list.import<-list(length=20)
for (i in 1:20)shapefiles.list.import[[i]]<-readOGR("buffer_shps/",shapefile.list[i])

### Load all upi shapefiles
## create an interation vector with the names endings
## creat a vector with the complete names
shapefile.list.upi<-vector(length=20) # create the vector

# create a list with the shapefile names
for (i in shapefile.iteration.vector)shapefile.list.upi[i/5]<-paste("_upi_buffer",i,sep="") # paste beginning and end of names
shapefile.list.upi[20]<-("_upi_buffer0") # _uus_buffer0 was lost so set it automatically
shapefile.list.upi
## import all shapefiles to the list
shapefiles.list.import.upi<-list(length=20)
for (j in 1:20)shapefiles.list.import.upi[[j]]<-readOGR("buffer_shps/",shapefile.list.upi[j])

### Load all indigenous shapefiles
## create an interation vector with the names endings
## creat a vector with the complete names
shapefile.list.ind<-vector(length=20) # create the vector

# create a list with the shapefile names
for (i in shapefile.iteration.vector)shapefile.list.ind[i/5]<-paste("_ind_buffer",i,sep="") # paste beginning and end of names
shapefile.list.ind[20]<-("_ind_buffer0") # _uus_buffer0 was lost so set it automatically
shapefile.list.ind

## import all shapefiles to the list
shapefiles.list.import.ind<-list(length=20)
for (k in 1:20)shapefiles.list.import.ind[[k]]<-readOGR("buffer_shps/",shapefile.list.ind[k])

### rearrange the shapefiles lists to have a better naming
shapefiles.list.import.2<-list()
shapefiles.list.import.ind.2<-list()
shapefiles.list.import.upi.2<-list()
shapefiles.list.import.2<-rev(shapefiles.list.import[1:19])
shapefiles.list.import.2[[20]]<-shapefiles.list.import[[20]]
shapefiles.list.import.ind.2<-rev(shapefiles.list.import.ind[1:19])
shapefiles.list.import.ind.2[[20]]<-shapefiles.list.import.ind[[20]]
shapefiles.list.import.upi.2<-rev(shapefiles.list.import.upi[1:19])
shapefiles.list.import.upi.2[[20]]<-shapefiles.list.import.upi[[20]]

# remove old shapefilelists
rm(shapefiles.list.import,shapefiles.list.import.ind,shapefiles.list.import.upi)

# find out the first year of the different type of protected areas
min(shapefiles.list.import.2[[20]]$DATE_CREAT,na.rm=T) # uus = 1927
min(shapefiles.list.import.ind.2[[20]]$DATE_CREAT,na.rm=T) # ind = 1926
min(shapefiles.list.import.upi.2[[20]]$DATE_CREAT,na.rm=T) # upi = 1896

### ----------C: construct parameters for the lineplot----------

# calculate sqkm for ucs cumulatice
for (i in 2:5)
  ucs_cum[paste((names(ucs_cum)[i]),"_sqkm",sep="")]<-ucs_cum[,i]/1000000

tail(ucs_cum)

# calculate the total protected area
ucs_cum$total_protected_area_sqkm<-ucs_cum$uus_m2_sqkm+ucs_cum$upi_m2_sqkm+ucs_cum$ind_m2_sqkm

# calculate uus+ind
ucs_cum$total_uus_ind_sqkm<-ucs_cum$uus_m2_sqkm+ucs_cum$ind_m2_sqkm

### --------- D: define viewports for the plots ----------
# define first plot area
vp1 <- viewport(x = 1.05, y = 0.133,
                height = 0.72, width = 0.72,
                just = c("right", "bottom"),
                name = "map"
)
# define second plot area
vp2 <- viewport(x = 0.05, y = 0.505,
                height = 0.3, width = 0.361,
                just = c("left", "bottom"),
                name="picture"
)
# define third plot area
vp3 <- viewport(x = 0.05, y = 0.12,
                height = 0.36, width = 0.361,
                just = c("left", "bottom"),
                name="lineplot"
)

# define fourth plot area
vp4 <- viewport(x = 0.5, y = 0.95,
                height = 0.1, width = 0.8,
                just = c("center", "center"),
                name="titel"
)

### ---------- E: define color schemes for the plots ----------
color.uus<-"#99FF33"
color.upi<-"#006600"
color.ind<-"#CC6600"

### --------- F: define the different themes for the plots ----------
## for additional information on how to manage them themes see the ggplot2 help
# http://docs.ggplot2.org/current/theme.html

# plot 1 = maps
theme.plot.1<-list(theme(panel.background = element_rect(fill = "lightskyblue2", color="black"),
                          panel.grid.minor = element_blank(),
                          panel.grid.major = element_blank(),
                          axis.text.x = element_text(size=25,colour="white"),
                          axis.text.y = element_text(size=25,colour="white"),
                          axis.title.x= element_text(size=25,colour="white"),
                          axis.title.y= element_text(size=25,colour="white"),
                          plot.title = element_text(size=75,vjust=3,family="DejaVu Sans Condensed",colour="white"),
                          plot.background = element_rect(fill="black",colour="black")
                         )
                   )
# plot 2 = lineplot
theme.plot.2<-list(theme(panel.grid.minor = element_blank(),
                         panel.grid.major = element_blank(),
                         plot.title = element_text(size=60,vjust=1,family="DejaVu Sans Condensed",colour="white"),
                         axis.text.x=element_text(size=25,angle=45,vjust=0.7,hjust=0.9,colour="white"),
                         axis.text.y=element_blank(),
                         axis.title.x=element_blank(),
                         axis.title.y=element_blank(),
                         panel.background=element_rect(fill = "black"),
                         legend.position=c(0,1),
                         legend.justification = c(0, 1),
                         legend.title=element_blank(),
                         legend.text=element_text(size=25,family="DejaVu Sans Condensed",colour="white"),
                         legend.key.height=unit(3,"line"),
                         legend.key.width=unit(3,"line"),
                         legend.background = element_rect(fill="black",colour="black"),
                         axis.ticks.length = unit(0.3, "cm"),
                         axis.line=element_line(),
                         axis.line.y = element_blank(),
                         plot.background = element_rect(fill="black",colour="black")
                         )
                   )

### --------- G: define additional parameters for the plots ----------
# the title
plot.title<-textGrob("Protected Areas in Brazil",
                     gp = gpar(fontfamily="DejaVu Sans Condensed",
                               fontsize = 90,
                               col="white"))

#### ---------- H: Construct the plotfunction ----------

plot.function.buffer<-function(my.path,
                               my.year,
                               data.unbuffered.1,
                               data.unbuffered.2,
                               data.unbuffered.3,
                               my.data.1,
                               my.data.2,
                               my.data.3,
                               my.buffer.size,
                               my.alpha,
                               my.picture)
{
  ### 1. fortify all shapefiles for the plot

  ## 1.1 fortify all existing protected areas until the year (unbufferd)
  # layer 1 unbuffered
  if(min(data.unbuffered.1@data$DATE_CREAT,na.rm=T)<=my.year+10000){
    tmp.data.unbuffered.1<-fortify(data.unbuffered.1[which
                                                     (data.unbuffered.1@data$DATE_CREAT<=my.year+10000),])}
  else{tmp.data.unbuffered.1<-"no area for first unbuffered layer"}
  # layer 2 unbuffered
  if(min(data.unbuffered.2@data$DATE_CREAT,na.rm=T)<=my.year+10000){
    tmp.data.unbuffered.2<-fortify(data.unbuffered.2[which
                                                     (data.unbuffered.2@data$DATE_CREAT<=my.year+10000),])}
  else{tmp.data.unbuffered.2<-"no area for second unbuffered layer"}
  # layer 3 unbuffered
  if(min(data.unbuffered.3@data$DATE_CREAT,na.rm=T)<=my.year+10000){
    tmp.data.unbuffered.3<-fortify(data.unbuffered.3[which
                                                     (data.unbuffered.3@data$DATE_CREAT<=my.year+10000),])}
  else{tmp.data.unbuffered.3<-"no area for third unbuffered layer"}

  ## 1.2 fortify protected areas for the respective year if there was a creation/data is available (with buffer)
  # layer 1 buffered
  if(nrow
     (my.data.1[[my.buffer.size]]@data[which
                                       (my.data.1[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                          my.data.1[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])>0
  )
  {tmp.data1<-fortify(my.data.1[[my.buffer.size]][which
                                                  (my.data.1[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                                     my.data.1[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])}else
                                                     {tmp.data1<-"no area for first layer"}
  # layer 2 buffered
  if(nrow
     (my.data.2[[my.buffer.size]]@data[which
                                       (my.data.2[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                          my.data.2[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])>0
  )
  {tmp.data2<-fortify(my.data.2[[my.buffer.size]][which
                                                  (my.data.2[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                                     my.data.2[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])}else
                                                     {tmp.data2<-"no area for second layer"}

  # layer 3 buffered
  if(nrow
     (my.data.3[[my.buffer.size]]@data[which
                                       (my.data.3[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                          my.data.3[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])>0
  )
  {tmp.data3<-fortify(my.data.3[[my.buffer.size]][which
                                                  (my.data.3[[my.buffer.size]]@data$DATE_CREAT>=my.year &
                                                     my.data.3[[my.buffer.size]]@data$DATE_CREAT<=my.year+10000),])}else
                                                     {tmp.data3<-"no area for third layer"}

  ### 2. construct the mapplot
  ## 2.1 add the base layers
  tmp.plot<-ggplot(,aes(long,lat))+
    geom_map(aes_string(map_id="id"),
             map=world.countries.fortified,
             data=world.countries.fortified,
             fill="grey80",
             colour="black",
             size=0.3)+
    geom_map(aes_string(map_id="id"),
             map=brasil.fortified,
             data=brasil.fortified,
             fill="lemonchiffon2",
             colour="black",
             size=0.3)

  ## 2.2 add all existing protected areas until the year unbufferd
  if (is.data.frame
      (tmp.data.unbuffered.1)==TRUE)
  {tmp.plot<-
     tmp.plot+geom_map(aes_string(map_id='id'),
                       map=tmp.data.unbuffered.1,
                       data=tmp.data.unbuffered.1,
                       fill=color.ind,
                       colour="black",
                       size=0.3)}
  # else{print(tmp.data.unbuffered.1)} # activate this option if you want to trace the layers within each map
  if (is.data.frame
      (tmp.data.unbuffered.2)==TRUE)
  {tmp.plot<-tmp.plot+geom_map(aes_string(map_id='id'),
                               map=tmp.data.unbuffered.2,
                               data=tmp.data.unbuffered.2,
                               fill=color.uus,
                               colour="black",
                               size=0.3)}
  # else{print(tmp.data.unbuffered.2)}
  if (is.data.frame
      (tmp.data.unbuffered.3)==TRUE)
  {tmp.plot<-tmp.plot+geom_map(aes_string(map_id='id'),
                               map=tmp.data.unbuffered.3,
                               data=tmp.data.unbuffered.3,
                               fill=color.upi,
                               colour="black",
                               size=0.3)}

  ## 2.3 add protected areas with buffer if data is available for the respective year
  # layer ind
  if (is.data.frame(tmp.data1)==TRUE){tmp.plot<-
                                        tmp.plot+geom_map(aes_string(map_id='id'),
                                                          map=tmp.data1, data=tmp.data1,
                                                          fill=color.ind,
                                                          colour="black",
                                                          size=0.3,
                                                          alpha=my.alpha)}
  # layer uus
  if (is.data.frame(tmp.data2)==TRUE){tmp.plot<-
                                        tmp.plot+geom_map(aes_string(map_id='id'),
                                                          map=tmp.data2, data=tmp.data2,
                                                          fill=color.uus,
                                                          colour="black",
                                                          size=0.3,
                                                          alpha=my.alpha)}
  # layer upi
  if (is.data.frame(tmp.data3)==TRUE){tmp.plot<-
                                        tmp.plot+geom_map(aes_string(map_id='id'),
                                                          map=tmp.data3, data=tmp.data3,
                                                          fill=color.upi,
                                                          colour="black",
                                                          size=0.3,
                                                          alpha=my.alpha)}

  ## 2.4 add limits, title, projection and theme configurations to the plot
  tmp.plot<-tmp.plot+
    xlim(-80,-30)+
    ylim(-33,5)+
    ggtitle(my.year/10000)+
    coord_map(projection="mercator")+
    theme.plot.1

  ### 3. construct the lineplot
  ## 3.1 plot
  tmp.plot2<-
    ggplot(data=ucs_cum[which(ucs_cum$year<my.year/10000),])+
    geom_ribbon(aes(x=year,# is going to be area for upi
                    ymin=total_uus_ind_sqkm,
                    ymax=total_protected_area_sqkm,
                    fill=color.upi,
                    color=color.upi))+
    geom_ribbon(aes(x=year,# is going to be area for uus
                    ymin=ind_m2_sqkm,
                    ymax=total_uus_ind_sqkm,
                    fill=color.uus,
                    color=color.uus))+
    geom_ribbon(aes(x=year,#is going to be area for ind
                    ymin=0,
                    ymax=ind_m2_sqkm,
                    fill=color.ind,
                    color=color.ind))+

    ## 3.2. area count as title
    ggtitle(
      paste(
        format(
          round(
            ucs_cum$total_protected_area_sqkm[which(iteration.string.plot==my.year)],
            digits=0), # together with round cuts the digits off
          big.mark=","), # together with format creates the thousend seperators)
        "km²"
      ))+
    ## 3.3 modify the plot
    # scales and colours
    scale_fill_manual(values=c(color.upi,color.uus,color.ind),
                      labels=c("  Strictly Protected Areas","  Multiple Use Reserves","  Indigenous Areas"))+
    scale_colour_manual(values=c(color.upi,color.uus,color.ind),
                        labels=element_blank(),
                        guide="none")+ # to surpress the legend of the outlines
    scale_x_continuous(breaks=seq(1920, 2010, by = 10),
                       limits=c(1920,2013), # set value limits of axis
                       expand=c(0, 0))+ # set visible limits of axis
    scale_y_continuous(breaks=NULL,
                       limits=c(0,2800000))+
    # add the themes
    theme.plot.2

  ### 4. Construct the plot with extra information
  tmp.plot3<-map.picture3.rastered

  ###  Combine the plots and make png outputs
  ## new png
  png(filename=paste(my.path,"map_",my.year,"_buffer",my.buffer.size,".png",sep=""), width=2560, height=1600)
  ## mapplot

  # mapplot
  grid.newpage()# clear plot area
  grid.rect(gp = gpar(fill = "black")) #show the plotting region (viewport extent)
  pushViewport(vp1) # enter vp1
  print(tmp.plot, newpage = FALSE) # plot1
  upViewport(1) # leave vp1

  # lineplot
  pushViewport(vp3)
  print(tmp.plot2, newpage = FALSE)
  upViewport(1)

  # mapinfo
  pushViewport(vp2)
  print(grid.draw(my.picture), newpage=FALSE)
  upViewport(1)

  # title
  pushViewport(vp4)
  print(grid.draw(plot.title), newpage=FALSE)
  upViewport(1)

  ## Turn the device off to finish the print
  dev.off()

  ### Finishing the function
  print(paste("Map completed for",my.year/10000,"with buffersize",my.buffer.size,"and opacity set to",my.alpha))
}

#### --------- I: Plot the data ----------

# define the path for the plot outputs
path.plot<-"r-outputs/"

# create a string to iterate over the input data ucs year by year
iteration.string.plot<-vector(length=length(1896:2013))
iteration.string.plot[1]<-18960000
for (i in 2:length(1896:2013))iteration.string.plot[i]<-iteration.string.plot[i-1]+10000
iteration.string.plot

# create a vector that contains the transparency/alpha values for each buffer
alpha.vector.one<-seq(from=0.5,to=1,length.out=20)

# plot all years and buffers (over 2gb size!)
for (j in 2:length(iteration.string.plot)){
  for (i in 19)
  {plot.function.buffer(path.plot,
                        iteration.string.plot[j],
                        shapefiles.list.import.ind.2[[20]],
                        shapefiles.list.import.2[[20]],
                        shapefiles.list.import.upi.2[[20]],
                        shapefiles.list.import.ind.2,
                        shapefiles.list.import.2,
                        shapefiles.list.import.upi.2,
                        i,
                        alpha.vector.one[i],
                        map.picture3.rastered)
  }
  Sys.sleep(3)
}

# ---------- J: Close session----------
save.image("r-workspaces/protected_areas_plotting_1.Rdata")

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s