R Notebook

This tutorial is to show how to produce a distribution map of poverty in Vietnam

Dataset is used for producing the map can be downloaded from this link

Point Data

Display spatial Point data using leaflet package. It is alo noted that all base maps can be found here

library(leaflet)
# Creating a random number of coordinates as below (data points are stored in a dataframe called df)
set.seed(123)

df<-data.frame(long=20,lat=runif(20))

leaflet(df) %>% addTiles()%>% addProviderTiles(providers$Esri.WorldImagery)%>% addCircles() 

# Alternative

leaflet(df) %>% addTiles() %>% addProviderTiles(providers$Thunderforest.SpinalMap) %>% addCircles(lng=~long,lat=~lat) 

*We can set a specific view for the map

# Display Ha Noi Vietnam

leaflet() %>% addTiles() %>% addProviderTiles(providers$OpenTopoMap) %>% setView(lng=105.834166,lat=21.027116,zoom=10)%>% addMarkers(lng=105.834166,lat=21.027116,popup="Welcome to Ha Noi Vietnam") # This is only fixed for initial display, while the user can zoom in and zoom out as desire

# Alternative

leaflet(options = leafletOptions(minZoom = 0,maxZoom = 10)) %>% addTiles() %>% addProviderTiles(providers$Esri.NatGeoWorldMap) %>% addMarkers(lng=105.834166,lat=21.027116,popup="Welcome to Ha Noi Vietnam", label="Welcome to Ha Noi Vietnam") # This will fix the viewer as you can't zoom in bigger than its current view

*Combining different tiles

leaflet() %>% addTiles()%>% addProviderTiles(providers$MtbMap) %>%
  addProviderTiles(providers$Stamen.TonerLines,
    options = providerTileOptions(opacity = 0.35)) %>%
  addProviderTiles(providers$Stamen.TonerLabels) %>% addMarkers(lng=105.834166,lat=21.027116)

*Makers Points

NZ<-read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/My_data.csv")

levels(NZ$ser) # Check the levels of Ser variable

## [1] "B" "C" "H" "P" "T"

head(NZ)

##   FID Field1       lat      lon ser
## 1   0      1 -37.74922 175.6729   T
## 2   1      2 -43.26288 170.3655   H
## 3   2      3 -37.67690 175.4695   P
## 4   3      4 -37.27856 175.4390   B
## 5   4      5 -37.25098 175.5982   C
## 6   5      6 -37.64233 175.3597   T

# The distribution of New Zealand species (unknown). This dataset belongs to Yunia (PhD student, Massey University)
leaflet(NZ) %>% addTiles() %>% addProviderTiles(providers$Thunderforest.Landscape) %>% addMarkers(lng=~lon,lat=~lat,popup = ~as.character(ser),label = ~as.character(ser))

*Makers Clusters

leaflet(NZ) %>% addTiles() %>% addMarkers(clusterOptions = markerClusterOptions()) 

## Assuming 'lon' and 'lat' are longitude and latitude, respectively

# Alternative using addCircles()

leaflet(NZ) %>% addTiles() %>% addCircleMarkers(lng=~lon,lat=~lat,radius=~c(2,3,4,5,6),color=~colorFactor(c("B", "C", "H", "P", "T"),domain=c("red","red","green","blue","yellow")),stroke=FALSE)

*AddCircles markers

pal <- colorNumeric("viridis", NULL)

leaflet(NZ) %>% addTiles() %>% addCircles(lng=~lon,lat=~lat,radius=~sqrt(Field1)*1000,weight=1,popup=~as.character(ser),label=~as.character(ser),color=~pal(Field1))

Getting online map

library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.2-8, (SVN revision 663)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/Tuyen/Documents/R/win-library/3.4/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/Tuyen/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5

vn<-getData(name="GADM",country="Vietnam",level=1)
# Due to some errors (dupplicates) in VN dataset, we will correct it. 
# Preprocess vn dataset so that they are consistent to each other

vn@data[vn@data$VARNAME_1=="Thua Thien - Hue","VARNAME_1"]<-"Hue"


vn@data[vn@data$VARNAME_1=="Ho Chi Minh City|Ho Chi Minh","VARNAME_1"]<-"Ho Chi Minh"

vn@data[vn@data$VARNAME_1=="Ha Noi City|Hanoi","VARNAME_1"]<-"Ha Noi"

vn@data[vn@data$VARNAME_1=="Hai Phong City|Haiphong","VARNAME_1"]<-"Hai Phong"

vn@data[vn@data$VARNAME_1=="Bac Kan|Bac Can","VARNAME_1"]<-"Bac Kan"

vn@data[vn@data$VARNAME_1=="Ba Ria - VTau|Ba Ria-Vung Tau","VARNAME_1"]<-"Ba Ria Vung Tau"

vn@data[vn@data$VARNAME_1=="Da Nang City|Da Nang","VARNAME_1"]<-"Da Nang"

vn@data[vn@data$VARNAME_1=="Dak Lak|Dac Lac","VARNAME_1"]<-"Dak Lak"
# There are two duplicates values and remove it as below

vn<-vn[-33,]

vn<-vn[-42,]

df<-read.csv("https://raw.githubusercontent.com/tuyenhavan/Statistics/Dataset/Poverty_Vietnam_2015.csv",sep=";") # This dataset was downloaded from the above-mentioned source.
df1<-df[order(df$Province),]

head(df1)

##          Province        Poor Poverty_Rate
## 1        An Giang  128,893.0          8.44
## 2 Ba Ria Vung Tau   95,205.0          1.91
## 3       Bac Giang   92,754.0         13.93
## 4         Bac Kan   81,592.0         29.40
## 5        Bac Lieu   74,313.0         15.55
## 6        Bac Ninh   65,374.0          3.53

vn1<-vn[order(vn$VARNAME_1),]


vn1$Province<-df$Province

vn1$Poverty_rate<-df$Poverty_Rate

vn1$Poor<-df$Poor

# Ready to create a map of poverty rate by geographic locations

pal<-colorNumeric("viridis", NULL)

leaflet() %>% addProviderTiles(providers$CartoDB.DarkMatter) %>% addPolygons(data=vn1,fillColor = ~pal(Poverty_rate),fillOpacity = 0.9,color="black",weight=0.5,popup=~as.character(Poverty_rate),label=~paste0(Province," Province",": ",Poverty_rate, " %")) %>% addLegend("bottomright",opacity=1,pal=pal,values=vn1$Poverty_rate,title="Poverty Rate by Province")