Mar 16, 2017

Mapping study area using R

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 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
#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.

8 comments:

  1. This blog aware me about different programs which can become very useful for our friends and kids. Few websites provide combined courses and few of the are separately for single subject. Glad to get this information.
    เรียน igcse math 

    ReplyDelete
  2. Share great information about your blog , Blog really helpful for us . We read your blog , share most useful information in blog . Thanks for share your blog here .
    สมัคร สอบ ged 

    ReplyDelete
  3. It is really a helpful blog to find some different source to add my knowledge. I came into aware of new professional blog and I am impressed with suggestions of author.
    sat สอบ

    ReplyDelete
  4. This is the first time I came to this blog and I found some relevant stuff here. Basically I keen to know new parameters of writing every-time and sometime it become really very hard to find such kind of platform.
    ติว sat

    ReplyDelete
  5. Hi there! I study your content from a long time. delhi top schools list I still prefer doing it because I nonetheless find out something hot and special. I am actually thankful that you submit that content.

    ReplyDelete

  6. I was reading some of your content on this website and I conceive this internet site is really informative shriji baba saraswati vidya mandir! Keep on putting up.

    ReplyDelete
  7. So it could be great to study subjects that you appreciate examining. An intriguing actuality is that life is the thing that we make it, so it is critical to pick the things that we truly need from life. study agent.

    ReplyDelete
  8. This article was written by a real thinking writer without a doubt.I agree many of the with the solid points made by the writer.I’ll be back day in and day for further new updates. điều kiện du học Nhật Bản

    ReplyDelete