Mapping Taiwan
Ed Liu
Preface
The purpose of this article is to simply teach you how to draw the basic map of Taiwan.
Cordially
Ed Liu
System language
If you use the MAC OS system, you may need the code below to set your local language for importing the chinese data
Sys.setlocale(category = "LC_ALL", locale = "zh_TW.big5")Packages we need
require(tidyverse)
require(sp)
require(maptools)
require(rgeos)
require(rgdal)
require(plyr)
require(broom)
require(ggrepel)Importing the sp object
The shapefiles in this tutorial came from Taiwan open data platform. Google something like “直轄市、縣市界線(TWD97經緯度)”.
The structure of spatial data: We use the readOGR() function, readOGR{rgdal} which can read OGR vector maps into spatial objects to import our mapping data due to the specifically spatial data structure. All the maps we plot by using the multiple points. Therefore, R would link every point in their order to draw the spatial polygon that represents the geographic region.
Firstly, we download the shp data from DATA.GOV.TW, put the file directory as " dsn="file route"". The layer should be TOWN_MOI_(dowmload date) or COUNTY_MOI_(download date) depending on which statistical area you want to plot (also the village level data follow the same way).
TWshp.ct <- readOGR(dsn = "/Users/liupochen_macbook/Desktop/mapdata201906180930",layer="COUNTY_MOI_1080617",
stringsAsFactors = F)
ggTWshp.ct<-tidy(TWshp.ct)
tmp <- data.frame(id=as.character(c(0:(dim(TWshp.ct@data)[1]-1))),
ct_code=TWshp.ct@data$COUNTYCODE,
county.eng=TWshp.ct@data$COUNTYENG,
stringsAsFactors = F)
ggTWshp.ct<-join(ggTWshp.ct,tmp,by="id")
TWshp <- readOGR(dsn = "/Users/liupochen_macbook/Desktop/mapdata201906180922",layer="TOWN_MOI_1080617",
stringsAsFactors = F)
ggTWshp.tn<-tidy(TWshp)
tmp2<-data.frame(id=as.character(c(0:(dim(TWshp@data)[1]-1))),
tn_code=TWshp@data$TOWNCODE,
ct_code=TWshp@data$COUNTYCODE,
town.eng=TWshp@data$TOWNENG,
stringsAsFactors = F)
tmp2<-(join(tmp2,tmp[,c(2,3)],by="ct_code"))
ggTWshp.tn<-join(ggTWshp.tn,tmp2,by="id")
rm(TWshp,tmp,tmp2)Mapping
Now, you would get the basic Taiwan spatial objects data in the data frame, ggTWshp.tn and ggTWshp.ct Then, let’s plot our world!
Country-level map
Please notice that if we want to plot not only the county level but also the district level, there should be given two shp file, COUNTY_MOI and TOWN_MOI into the graphic grammar, ggplot2
ggplot()+
geom_polygon(aes(x = long, y=lat, group=group),data=ggTWshp.ct,
color="#d0e8f2",fill="#4fd67e", lwd=.3)+
geom_polygon(aes(x = long, y=lat, group=group),data=ggTWshp.tn,
color="#d0e8f2",alpha=.45, lwd=.07)+
coord_map(xlim=c(119.7,122.2),ylim=c(21.7,25.4))
For visualizability, coord_map() can help us to fix the coordinate system.
Most of the time, we want to zoom in for presenting the district level, even the township level. Take “Kaohsiung City” for example, we can use dplyr::filter() to filter out the district from the ggTWshp then plot the high-resolution map for Kaohsiung City. (Of the Kaohsoung map, I limit the coordinate system to focusing on the mainland of Taiwan by xlim for longitude and ylim for latitude.)
ggKao.ct<- ggTWshp.ct%>%filter(county.eng=="Kaohsiung City")
ggKao.tn<- ggTWshp.tn%>%filter(county.eng=="Kaohsiung City")
#Mapping
ggplot()+
geom_polygon(aes(x = long, y=lat, group=group),data=ggKao.ct,
color="#d0e8f2",fill="#4fd67e", lwd=.3)+
geom_polygon(aes(x = long, y=lat, group=group),data=ggKao.tn,
color="#d0e8f2",alpha=.45, lwd=.07)+
coord_map(xlim=c(120.1,121.1),ylim=c(22.4,23.5)) 
Finally, for the basic mapping, if we want to plot the point on the map, we could simplly plus +geom_point() in our plotting grammar. But, hold on! Just quickly have a peek on our Kaohsiung map, can you simply point out every district names? Maybe for the majority of Taiwanese could not. So if we label each district name on the map, that would be helpful.
Let’s look deeper in it!
I am going to present the most simple way to calculate the specific area center to present how to add the geom_point() function. Here I calculate the center of all the points for the particular polygon by the function which I give it into cnames. Therefore, add the other layer on the original mapping grammar for labeling the district names. geom_label_repel help us to point out the center of each district in Kaohsung and link the district name box from the district center. (about ggrepel::geom_label_repel(), the package allows everyone to customizely design the decoration. Function detail would be described in “R Documentation”)
cnames <- aggregate(cbind(long, lat) ~ town.eng, data=ggKao.tn,
FUN=function(x)mean(range(x)))
ggplot()+
geom_polygon(aes(x = long, y=lat, group=group),data=ggKao.ct,
color="#d0e8f2",fill="#4fd67e", lwd=.3)+
geom_polygon(aes(x = long, y=lat, group=group),data=ggKao.tn,
color="#d0e8f2",alpha=.45, lwd=.07)+
geom_label_repel(data=cnames,aes(x=long,y=lat,label=town.eng),
size=1.32,
col="#0D859C",
box.padding = 0.55,
point.padding = 0.3,
label.size = 0.35,
segment.color = "#0D859C",
alpha=0.6,
position="identity"
)+
coord_map(xlim=c(120.1,121.1),ylim=c(22.4,23.5)) 
Choropleth map
When it comes to mapping, the choropleth map should not be left behind. First of all, we may need some prepare. The map theme below is which I usually use. (As you see below, you can also design for yourself!)
theme_bare <- theme(
axis.line = element_blank(),
panel.grid.major.x = element_line(color="#53dba5",size =0.1),
panel.grid.major.y = element_line(color="#53dba5",size=0.1),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
legend.text=element_text(size=7),
legend.title=element_text(size=8),
panel.background = element_blank(),
panel.border = element_rect(colour = "gray", fill=NA, size=0)
)The choropleth map can elaborate the area statistics in an intuitive way for looking at. However, how the choropleth map works? And how to plot the choropleth map? Each area in the choropleth map could represent their own statistical value. For example, the incidence proportion or odd ratio. Therefore, we fill the specific area by different colors to represent the scale of the statistical values. In the data manipulating process, we have to create a new column for the statistics that you want to plot on the map. Firstly, I create an example dataset, cnames_eg (from cnames), where the incidence risk representing the statistics for mapping. There leaves only one step should be done, merge the values into the spatial object (here refers as ggKao.ct or ggKao.tn, depends on which statistical areas you would like to communicate to your audience) then we can plot it by filling incidence risk in the particular polygons.
About merge
dplyr::left_join() will put the incidence risk values in every points for mapping the Kaohsiung districts by the column names town.eng and without changing the order of the map points.
cnames_eg<- cnames%>%mutate(incidence_risk=c(1:38))
ggKao.tn.cp<-dplyr::left_join(ggKao.tn,cnames_eg[,c(1,4)], by="town.eng")Fill the polygon with the statistics
For the last step, the only thing you have to do is assigning the fill= in geom_polygon(aes()) to the statistics, here incidence_risk.
#mapping
ggplot()+
geom_polygon(data=ggKao.tn.cp, aes(x=long, y=lat, group=group,fill=incidence_risk),size=1.5)+
coord_map(xlim=c(120.1,121.1),ylim=c(22.4,23.5))+
scale_fill_gradient(high = "#e34a33",
low = "#fee8c8",
guide = "colorbar",
name="Incidence risk",
limit=c(0,40))+
theme(legend.key.height= unit(2.5,"mm"))+
theme(legend.background=element_blank())+
theme(legend.key=element_blank())+
theme(legend.position = c(1,0),
legend.justification = c(1,0),
legend.direction = "horizontal",
legend.box = "horizontal")+
theme_bare+
labs(title="Choropleth map example")
There are many ways to create a map, the method of this article is one of the simpliest one. The last thing I want to mention is that the projection of the coordinates will effect how the map look like, however, if we want to calculate some spatial indicators which need the real distance, we should take attention at the point coordinates.
Happy mapping!
Ed
Supplement {Google Map Web Crawler}
require(rjson)
require(httr)
#### API key from https://developers.google.com/maps/documentation/geocoding/start
getLatLng <- function(address){
address <- as.character(address)
address <- iconv(address,"big5","utf8")
urlData <- GET(paste0("https://maps.googleapis.com/maps/api/geocode/json?language=zh-TW&address=", URLencode(address),"&key= Your google API Key "))
jsonResult <- rjson::fromJSON(rawToChar(urlData$content))
Sys.sleep(1)
if(jsonResult$status != "OK"){
print("Google geocode API Error")
return("error")
}
print("LatLng Got")
lat <<- jsonResult$results[[1]]$geometry$location$lat
lng <<- jsonResult$results[[1]]$geometry$location$lng
return(c(lat, lng))
}