To prepare a study area map in a professional way, or at publication standard, it takes a great deal of time, at least for me. I have been using GIS tools (such as QGIS and Grass GIS, Jgrasstools) to prepare my study area. But, Now I started to change to R. The advantages of using R is that we have more advanced graphical and map package that we can use, and once the R script are developed, it is easier to reproduce other maps for the next time. Here I used various packages to produce decent study area map that can be used for publication. I used many helps from various stackoverflow answers. Here I will shows you how to do this using a particular basin in Afghanistan.
First lets load all the necessary r packages
library(rgdal)
library(ggplot2)
library(reshape2)
library(maptools)
## Checking rgeos availability: TRUE
library(sp)
library(rgeos)
## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-2
## Polygon checking: TRUE
library(ggplot2)
library(raster)
library(gridExtra)
library(raster)
library(grid)
Country map
AF0<-getData("GADM", country="AF", level=0)
AF0.f <- fortify(AF0, region = "ISO")
The study basin
mybasin<-readOGR(dsn="/Users/administrator/Documents/Afganistan", layer="mybsinMask")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/administrator/Documents/Afganistan", layer: "mybsinMask"
## with 1 features
## It has 6 fields
utms <- CRS("+init=epsg:32642")
proj4string(mybasin)<-utms
mybasin2 <- spTransform(mybasin, CRS( "+init=epsg:4326"))
mybasin2.f <- fortify(mybasin2, region = "area")
More information of the study basin
What I ma doing is hydrological anlayis, so subbasin partition is necessary. Here, I will plot the subbasin polygon maps as follows, as showsn for the basin polygon.
polysub<-readOGR(dsn="/Users/administrator/Documents/Afganistan", layer="Subbasin_final")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/administrator/Documents/Afganistan", layer: "Subbasin_final"
## with 15 features
## It has 7 fields
utms <- CRS("+init=epsg:32642")
proj4string(polysub)<-utms
polysub2 <- spTransform(polysub, CRS( "+init=epsg:4326"))
polysub2.f <- fortify(polysub2, region = "netnum")
Similarly we do the same for the river netwrok shapefile.
Network<-readOGR(dsn="/Users/administrator/Documents/Afganistan", layer="Network_final")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/administrator/Documents/Afganistan", layer: "Network_final"
## with 15 features
## It has 6 fields
utms <- CRS("+init=epsg:32642")
proj4string(Network)<-utms
Network2 <- spTransform(Network, CRS( "+init=epsg:4326"))
Network2.f <- fortify(Network2, region = "netnum")
In the following, we have the dem of the study area, and here we would like to show it on the map
Farah_raster<-raster("/Users/administrator/Documents/Afganistan/20170220070801_215846837.asc")
Farah_raster<-mask(Farah_raster, mybasin2)
Farah_raster2 <- rasterToPoints(Farah_raster)
Farah_raster3 <- data.frame(Farah_raster2)
colnames(Farah_raster3) <- c("long", "lat", "Elev")
Plotting using ggplot2
Now, we can use those data frames to use ggplot functionalities to plot high quality maps. First we plot the basin, subbasin polyogn, and the river network overlaid on the dem elevation for the study basin. Then, we will use the inset function that can be inserted in the figure to show the country location for those who are not familiar to the study area (basin).
p1<-ggplot()+ geom_raster(data=Farah_raster3, aes(long,lat, fill=Elev)) +
geom_polygon(data=polysub2, aes(long,lat, group=group), fill="grey40", colour="black",size=0.4, alpha=0)+
geom_path(data=Network2, aes(long,lat, group=group), colour="blue")+
coord_equal()+theme_bw()+xlab("")+ylab("") +
scale_fill_gradientn(name="Elevation (m)", colours = terrain.colors(10)) +
theme(legend.position = c(0.13, 0.75))
## Regions defined for each Polygons
The second plot i.e. the country shapefile.
p2<-ggplot()+geom_polygon(data=AF0, aes(long,lat,group=group),colour="grey10",fill="#fff7bc")+
geom_line(data=mybasin2, aes(long,lat, group=group))+
coord_equal()+theme_bw()+labs(x=NULL,y=NULL)+
theme(axis.text.x =element_blank(),axis.text.y= element_blank(), axis.ticks=element_blank(),axis.title.x =element_blank(),
axis.title.y= element_blank())
## Regions defined for each Polygons
## Regions defined for each Polygons
And finaly we can use the grid to combine the two plots and inset the country map inside the study area figure frame.
grid.newpage()
v1<-viewport(width = 1, height = 1, x = 0.5, y = 0.5)
v2<-viewport(width = 0.28, height = 0.3, x = 0.84, y = 0.38)
print(p1,vp=v1)
print(p2,vp=v2)
Finally we obtained good quality study area map.