How to create a map in R with ggplot2

Introduction

It is possibe to create maps in ggplot2 to ensure reproducibility. # Loading packages

The codes bellow will load the required packages used for the analyses. The function ipak is set to install the packages if they are not yet installed and then load them. To know which package was used for a particular function, simply press F2 on that function and the details will open in a new windows.

package_list=c('dendextend','ggdendro','ggplot2','ggpubr','ggthemes',
               'indicspecies','MuMIn','plyr','dplyr','png','questionr',
               'adehabitatHS', 'knitr','aod','car','cluster',
               'corrr','cowplot','DescTools','factoextra','fpc',
               'interactions','openxlsx','reshape2','rsq',
               'vegan','ggrepel','here','ggstar','pivottabler',
               'optpart','collapse', 'tmap','tmaptools','sf','ggmap',
               'rcartocolor', 'cowplot','grid','ggspatial','ggsn')

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

ipak(package_list)

##   dendextend     ggdendro      ggplot2       ggpubr     ggthemes indicspecies 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##        MuMIn         plyr        dplyr          png    questionr adehabitatHS 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##        knitr          aod          car      cluster        corrr      cowplot 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##    DescTools   factoextra          fpc interactions     openxlsx     reshape2 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##          rsq        vegan      ggrepel         here       ggstar  pivottabler 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##      optpart     collapse         tmap    tmaptools           sf        ggmap 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##  rcartocolor      cowplot         grid    ggspatial         ggsn 
##         TRUE         TRUE         TRUE         TRUE         TRUE

Creating the map of the study area

Africa_sf <- st_read("C:/Harvard_Dataverse/African_Continent.shp")

## Reading layer `African_Continent' from data source 
##   `C:\Harvard_Dataverse\African_Continent.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -25.36055 ymin: -34.822 xmax: 63.49576 ymax: 37.34041
## Geodetic CRS:  WGS 84

Botany_sf <- st_read("C:/Harvard_Dataverse/BotanicalPlots_All_Types.shp")

## Reading layer `BotanicalPlots_All_Types' from data source 
##   `C:\Harvard_Dataverse\BotanicalPlots_All_Types.shp' using driver `ESRI Shapefile'
## Simple feature collection with 575 features and 11 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 282080 ymin: 370236 xmax: 297261 ymax: 381750
## Projected CRS: WGS 84 / UTM zone 33N

Dja_sf <- st_read("C:/Harvard_Dataverse/Dja_river2.shp")

## Reading layer `Dja_river2' from data source `C:\Harvard_Dataverse\Dja_river2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 240334.8 ymin: 344314.4 xmax: 333468.3 ymax: 417074.9
## Projected CRS: WGS 84 / UTM zone 33N

Grid_sf <- st_read("C:/Harvard_Dataverse/Grid_1000x1000.shp")

## Reading layer `Grid_1000x1000' from data source 
##   `C:\Harvard_Dataverse\Grid_1000x1000.shp' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 282000 ymin: 370000 xmax: 298000 ymax: 382000
## Projected CRS: WGS 84 / UTM zone 33N

CMR_sf <- st_read("C:/Harvard_Dataverse/CMR_Outline_2011_dd.shp")

## Reading layer `CMR_Outline_2011_dd' from data source 
##   `C:\Harvard_Dataverse\CMR_Outline_2011_dd.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 8.491707 ymin: 1.654588 xmax: 16.1891 ymax: 13.08015
## Geodetic CRS:  WGS 84

## Extent of Cameroon in Africa
CMR_sf_32633 = st_transform(CMR_sf, crs = 32633)

ggm1 = ggplot() + 
  geom_sf(data = Africa_sf, fill = "white", color = "gray70") + 
  geom_sf(data = CMR_sf, fill = "gray70", color = "gray70") +
  theme_void()


## Extent of study area in Cameroon
Grid_sf_bb = st_as_sfc(st_bbox(Grid_sf))
Grid_sf_bb2 = st_buffer(Grid_sf_bb, dist = 1000)
ggm2 = ggplot() + 
  geom_sf(data = CMR_sf_32633, fill = "white", color = "gray70") + 
  geom_sf(data = Grid_sf_bb2, fill = NA, color = "blue", size = 1.1) +
  theme_void()


## Map of the study area
# the function 'coord_sf' is used to limit the extent of our main map to the 
# range of the frame in the inset map
Colors_used <- c("#E58606","#5D69B1","#52BCA3","#99C945","#24796C","#ED645A","#CC3A8E","#A5AA99")
Shapes_used <- c(7,9,6,1,24,23,22,21)
Botany_sf$Type2=ifelse(Botany_sf$PlotType=="Nesting plot","Sleeping plot", "Systematic plot")
Botany_sf$Variable=paste0(Botany_sf$Type2,", ",Botany_sf$Habitat11)

ggm3 = ggplot() + 
  geom_sf(data = Grid_sf, fill = "#F5F5DC") +
  geom_sf(data = Botany_sf, aes(shape=Variable, fill=Variable)) +
  theme_bw() +
  theme(legend.position = c(0.5, 0.9),
        legend.direction = "horizontal",
        legend.key.width = unit(10, "mm"),
        plot.background = element_rect(fill = NA),
        legend.background = element_rect(fill = NA)) +
  scale_shape_manual(values = Shapes_used) +
  scale_color_manual(values = Colors_used)+
  guides(shape=guide_legend(nrow=4, byrow=FALSE))+
  theme(legend.title=element_blank())+labs(x = NULL, y = NULL)+
  coord_sf(xlim = c(st_bbox(Grid_sf_bb2)[c(1)]-3500,
                    st_bbox(Grid_sf_bb2)[c(3)]),
           ylim = c(st_bbox(Grid_sf_bb2)[c(2)]-1500,
                    st_bbox(Grid_sf_bb2)[c(4)])+3000)+
  scalebar(Grid_sf, dist = 2, dist_unit = "km",
           transform = FALSE, model = "WGS84",
           location="topleft")+
  annotation_north_arrow(which_north = "true")

## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## i Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Joining the maps
Study_Area_map = ggdraw() +
  draw_plot(ggm3) +
  draw_plot(ggm1, x = 0.05, y = 0.65, width = 0.27, height = 0.27)+
  draw_plot(ggm2, x = 0.001, y = 0.3, width = 0.35, height = 0.35)




setwd("C:/Harvard_Dataverse/Results")
## Plots for clusters
tiff(file = "Study_Area_map.tif", width = 20, height = 15, units = "cm", pointsize = 10,
     bg = "white", res = 300)
Study_Area_map
dev.off()

## png 
##   2

HAb_dir=setwd("C:/Harvard_Dataverse/Results")

ggsave("Study_Area_map.pdf",Study_Area_map,  path = HAb_dir,
       scale = 1, width = 20, height = 15, units = "cm",
       dpi = 600)

ggsave("Study_Area_map.eps",Study_Area_map,  path = HAb_dir,
       scale = 1, width = 20, height = 15, units = "cm",
       dpi = 600)

Location of the study area. OSF, old secondary forest; YSF, young secondary forest; SW, swamp; RIP, riparian forest