Search code examples
rggplot2ggmapr-sfr-sp

sf object is not properly overlaid on ggmap layer in r


I am trying to plot sf object over ggmap terrain layer in R. I am using the following code

library(ggmap)
library(sf)
library(tidyverse)

#Downloading data from DIVA GIS website
get_india_map <- function(cong=113) {
  tmp_file <- tempfile()
  tmp_dir  <- tempdir()
  zp <- sprintf("http://biogeo.ucdavis.edu/data/diva/adm/IND_adm.zip",cong)
  download.file(zp, tmp_file)
  unzip(zipfile = tmp_file, exdir = tmp_dir)
  fpath <- paste(tmp_dir)
  st_read(fpath, layer = "IND_adm2")
}
ind <- get_india_map(114)

#To view the attributes & first 3 attribute values of the data
ind[1:3,]

#Selecting specific districts
Gujarat <- ind %>% 
  filter(NAME_1=="Gujarat") %>%
  mutate(DISTRICT = as.character(NAME_2)) %>%
  select(DISTRICT)

#Added data to plot
aci <- tibble(DISTRICT=Gujarat$DISTRICT,
       aci=c(0.15,0.11,0.17,0.12,0.14,0.14,0.19,0.23,0.12,0.22,
                         0.07,0.11,0.07,0.13,0.03,0.07,0.06,0.04,0.05,0.04,
                         0.03,0.01,0.06,0.05,0.1))

Gujarat <- Gujarat %>% left_join(aci, by="DISTRICT")

#Plotting terrain layer using ggmap
vt <- get_map("India", zoom = 5, maptype = "terrain", source = "google")
ggmap(vt)

#Overlaying 'sf' layer
ggmap(vt) + 
  geom_sf(data=Gujarat,aes(fill=`aci`), inherit.aes=F, alpha=0.9) + 
  scale_fill_distiller(palette = "Spectral")

which returns me

enter image description here

As you can see from the plot the sf layer is not overlaid properly on the ggmap terrain layer. How to properly overlay the sf layer on the ggmap terrain layer?

But When I am using sp object in place of sf object the polygon fits properly on ggmap like

library(sp)
# sf -> sp
Gujarat_sp <- as_Spatial(Gujarat) 

viet2<- fortify(Gujarat_sp)
ggmap(vt) + geom_polygon(aes(x=long, y=lat, group=group), 
                         size=.2, color='black', data=viet2, alpha=0) + 
  theme_map() + coord_map()

enter image description here

But I don't know how to fill the geom_polygon according to aci?


Solution

  • I have solved the problem using sp package after taking help from this, this and this. The problem was that after applying fortify on the spatial data, only polygon information was there without the attribute data. So, I have merged the attribute data to the fortified polygon. Here is the complete code

    library(ggmap)
    library(sf)
    library(tidyverse)
    library(sp)
    
    #Downloading data from DIVA GIS website
    get_india_map <- function(cong=113) {
      tmp_file <- tempfile()
      tmp_dir  <- tempdir()
      zp <- sprintf("http://biogeo.ucdavis.edu/data/diva/adm/IND_adm.zip",cong)
      download.file(zp, tmp_file)
      unzip(zipfile = tmp_file, exdir = tmp_dir)
      fpath <- paste(tmp_dir)
      st_read(fpath, layer = "IND_adm2")
    }
    ind <- get_india_map(114)
    
    #To view the attributes & first 3 attribute values of the data
    ind[1:3,]
    #To plot it
    plot(ind["NAME_2"])
    
    #Selecting specific districts
    Gujarat <- ind %>% 
      filter(NAME_1=="Gujarat") %>%
      mutate(DISTRICT = as.character(NAME_2)) %>%
      select(DISTRICT)
    
    #Creating some data
    aci <- tibble(DISTRICT=Gujarat$DISTRICT,
           aci=c(0.15,0.11,0.17,0.12,0.14,0.14,0.19,0.23,0.12,0.22,
                             0.07,0.11,0.07,0.13,0.03,0.07,0.06,0.04,0.05,0.04,
                             0.03,0.01,0.06,0.05,0.1))
    
    vt <- get_map("India", zoom = 5, maptype = "terrain", source = "google")
    ggmap(vt)
    
    # sf -> sp
    Gujarat_sp <- as_Spatial(Gujarat) 
    
    # fortify the shape file
    viet2<- fortify(Gujarat_sp, region = "DISTRICT")
    
    # merge data
    map.df <- left_join(viet2, aci, by=c('id'='DISTRICT')) 
    
    #Plotting
    ggmap(vt) + geom_polygon(aes(x=long, y=lat, group=group, fill = `aci`), 
                             size=.2, color='black', data=map.df, alpha=0.8) + 
      theme_map() + coord_map() + scale_fill_distiller(name = "ACI", palette = "Spectral")
    

    enter image description here

    To make discrete classes following code may be used

    #For discrete classes
    map.df$brks <- cut(map.df$aci, 
                       breaks=c(0, 0.05, 0.1, 0.15, 0.2, 0.25), 
                       labels=c("0 - 0.05", "0.05 - 0.10", "0.10 - 0.15", 
                                "0.15 - 0.20", "0.20 - 0.25"))
    
    # Mapping with the order of colour reversed 
    ggmap(vt) + geom_polygon(aes(x=long, y=lat, group=group, fill = brks), 
                         size=.2, color='black', data=map.df, alpha=0.8) +
      theme_map() + coord_map() + 
      scale_fill_brewer(name="ACI", palette = "Spectral", direction = -1)
    

    enter image description here