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
#we can use getData fiunction to access avaliable geographic data
AF0<-getData("GADM", country="AF", level=0) # download country level 0 map
## In case if we want to chnag it to data frame using fortify, reshape2 function
AF0.f <- fortify(AF0, region = "ISO")The study basin
#Now lets read the basin extracted in other GIS environment using rgdal package
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#Lets convert from the projection of the original map to the latlong
utms <- CRS("+init=epsg:32642")
proj4string(mybasin)<-utms
mybasin2 <- spTransform(mybasin, CRS( "+init=epsg:4326"))
#Now we need to merge this with other shapefile we wanted to plot, and for that we need to change shapefile to dataframe. For that we used fortify.
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 fieldsutms <- 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 fieldsutms <- 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
#First we reed the raster using the raster function
Farah_raster<-raster("/Users/administrator/Documents/Afganistan/20170220070801_215846837.asc")
#here we like to take only the dem for the basin, and delet map outside the basin
Farah_raster<-mask(Farah_raster, mybasin2)
#convert the raster to points for plotting
Farah_raster2 <- rasterToPoints(Farah_raster)
#Make the points a dataframe for ggplot so that we can merg with the polygon iformation
Farah_raster3 <- data.frame(Farah_raster2)
#Make appropriate column headings, wich is common with the other dataframe
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).
# here is the basin, subbasin, river network, and elevatin Map
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) #plot area for the main map
v2<-viewport(width = 0.28, height = 0.3, x = 0.84, y = 0.38) #plot area for the inset map
print(p1,vp=v1)
print(p2,vp=v2)
Finally we obtained good quality study area map.
