#################################################
# E04: R語言畫地圖 (Maps)                       #
# 吳漢銘 國立政治大學統計學系                   #
# https://hmwu.idv.tw                           #
#################################################

##############################################################
# Setup                                                      #
##############################################################

Sys.setlocale(category = "LC_ALL", locale = "zh_TW.big5")

##############################################################
# 每個人都要有 Google 帳號: 例如: myname@gmail.com           #
##############################################################
# step1: Get API Key
# https://cloud.google.com/maps-platform/
# https://developers.google.com/maps/documentation/maps-static/get-api-key
# Step2: Enable Billing on the Google Cloud Project at 
# https://console.cloud.google.com/project/_/billing/enable 


# 5/54
library(ggmap)
register_google(key = "AIzaSyCuYcvrytmKLGN", write = TRUE) 
tw.xy <- geocode("Taiwan")
tw.xy
has_google_key()
google_key()

tw.map <- get_map(location = 'Taiwan', zoom = 7, language = "zh-TW")


# 7/54
library(RgoogleMaps)
TaiwanMap <- GetMap(center = c(lat = 23.58, lon  = 120.58), zoom  = 7, destfile = "Taiwan1.png")
TaiwanMap <- GetMap(center = c(lat = 23.58, lon  = 120.58), zoom = 10, destfile = "Taiwan2.png", maptype = "terrain")


# 8/54
my.lat <- c(25.175339, 25.082288, 25.042185, 25.046254)
my.lon <- c(121.450003, 121.565481, 121.614548, 121.517532)
bb = qbbox(my.lat, my.lon)
print(bb)
MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "my.png", maptype = "roadmap")
My.markers <- cbind.data.frame(lat = my.lat, lon = my.lon)
tmp <-  PlotOnStaticMap(MyMap, lat = My.markers[,"lat"], lon = My.markers[,"lon"], destfile = "my.png", cex = 2.5, pch = 20, col = 1:4, add = F)


# 9/54
png("my2.png", 640, 640)
tmp <-  PlotOnStaticMap(MyMap, lat = My.markers[,"lat"], lon = My.markers[,"lon"], cex = 2.5, pch = 20, col = 1:4, add = F)
tmp <-  PlotOnStaticMap(MyMap, lat = My.markers[,"lat"], lon = My.markers[,"lon"], col = "blue", add = T, FUN = lines, lwd = 2)
dev.off()


# 10/54
my.lat <- c(25.175339, 25.14362, 24.942605)
my.lon <- c(121.450003, 121.501768, 121.368381)
bb = qbbox(my.lat, my.lon)
print(bb)
MyMap <- GetMap.bbox(bb$lonR, bb$latR, destfile = "my.png", maptype = "roadmap")

My.markers <- cbind.data.frame(lat = my.lat, lon = my.lon)
tmp <- PlotOnStaticMap(MyMap, lat = My.markers[,"lat"], lon = My.markers[,"lon"], 
                       destfile = "my.png", cex = 2.5, pch = 18:10, col = 1:3, add = F)
TextOnStaticMap(MyMap, lat = My.markers[,"lat"]+0.01, lon = My.markers[,"lon"], 
                labels = c("我家", "復興高中", "國立臺北大學三峽校區"),  add = T)


# install EBImage
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("EBImage")

library(EBImage)
ntpu <- readImage("data/NTPUcolorlogo.jpg")
loc <- LatLon2XY.centered(MyMap, lat = My.markers[3, 1], lon = My.markers[3, 2])
rasterImage(ntpu, loc[[1]], loc[[2]]+30, loc[[1]]+50, loc[[2]]+80)
Fuxing <- readImage("data/Fuxinglogo.jpg")
loc <- LatLon2XY.centered(MyMap, lat = My.markers[2, 1], lon = My.markers[2, 2])
rasterImage(Fuxing, loc[[1]], loc[[2]]+30, loc[[1]]+50, loc[[2]]+80)



# 11/54
library(maps); library(maptools); library(mapdata); library(mapproj)
layout(matrix(c(1,1,1,0,2,0), ncol = 2), widths = c(10, 1), heights = c(1, 10, 1))
map("world2Hires", xlim = c(118, 123), ylim = c(21, 26))
data <- read.table("20140714-weather.txt", sep = "\t", header = TRUE, row.names = 1)
x <- data$TEMP
tm <- floor((100-1)/(max(x)-min(x))*(x-min(x)) + 1)
used.col <- heat.colors(100)[tm]
points(data$lon, data$lat, pch = 15, col = used.col)
text(data$lon, data$lat, labels = row.names(data))
title("20140714, 晚上8時各地溫度")
par(mar = c(1,1,1,1))
image(t(matrix(c(1:100), ncol = 1)), 
      col = heat.colors(100), xaxt = "n", yaxt = "n")
axis(LEFT <- 2, at = tm/100, 
     labels = as.character(x), cex.axis = 1)


# 14/54
tw.map <- get_map(location = 'Taiwan', zoom = 7)
ggmap(tw.map)

tw.map.zh <- get_map(location = 'Taiwan', zoom = 7,
                     language = "zh-TW")
ggmap(tw.map.zh)


# 15/54
tw.map.ntpu <- get_map(location = c(lon = 121.374925, lat = 24.943403),
                       zoom = 10, language = "zh-TW", maptype = "terrain")
ggmap(tw.map.ntpu)

map <- get_map(location = c(lon = 121.374925, lat = 24.943403),
               zoom = 10, language = "zh-TW", maptype = "toner-lite")
ggmap(map)



# 16/54
uv <- read.csv("data/UV_20191016130309.csv")
head(uv)

# 經緯度(度分秒)  = > 度數
lon.deg <- sapply((strsplit(as.character(uv$WGS84Lon), ",")), as.numeric)
lon.deg
uv$lon <- lon.deg[1, ] + lon.deg[2, ]/60 + lon.deg[3, ]/3600
lat.deg <- sapply((strsplit(as.character(uv$WGS84Lat), ",")), as.numeric)
uv$lat <- lat.deg[1, ] + lat.deg[2, ]/60 + lat.deg[3, ]/3600

tw.map <- get_map(location = 'Taiwan', zoom = 7)
ggmap(tw.map) + 
  geom_point(data = uv, aes(x = lon, y = lat, size = UVI), color = "purple")


# 17/54
ggmap(tw.map) + 
  geom_point(data = uv, aes(x = lon, y = lat, size = UVI, color = UVI)) +
  scale_color_continuous(low = "yellow", high = "red") +
  facet_grid(~PublishAgency) +
  guides(size = FALSE)


# 18/54
tpe.map.zh11 <- get_map(location = 'Taipei', zoom = 11, maptype = "roadmap", language = "zh-TW")
ggmap(tpe.map.zh11)
tpe.map.zh12 <- get_map(location = 'Taipei', zoom = 12, maptype = "roadmap", language = "zh-TW")
ggmap(tpe.map.zh12)
tpe.map.zh15 <- get_map(location = 'Taipei', zoom = 15, maptype = "roadmap", language = "zh-TW")
ggmap(tpe.map.zh15)
tpe.map.zh21 <- get_map(location = 'Taipei', zoom = 21, maptype = "roadmap", language = "zh-TW")
ggmap(tpe.map.zh21)


# 19/54
head(crime)
dim(crime)
summary(crime$offense)
library(dplyr)
rapes <- filter(crime, offense == "rape") %>%
  select(date, offense, address, lon, lat)
head(rapes)
dim(rapes)


# 20/54
houston_center <- geocode("Houston, TX")

# register_google(key = "AIzaSyCuYcvrytmKLGNxxx", write = TRUE) 
houston_center <- geocode("Houston, TX")
houston_center

has_google_key()
google_key()
ggmap_show_api_key()

houston_map <- get_map(houston_center, zoom = 13, maptype = "roadmap")
ggmap(houston_map)


# 21/54
ggmap(houston_map,
      base_layer = ggplot(data = rapes, aes(x = lon, y = lat))) +
  geom_point(color = "red", size = 3, alpha = 0.5)

ggmap(houston_map) +
  geom_point(data = rapes, aes(x = lon, y = lat), 
             color = "red", size = 3, alpha = 0.5)

ggmap(houston_map) +
  geom_point(data = rapes, aes(x = lon, y = lat), 
             color = "red", size = 3, alpha = 0.5) +
  geom_density2d(size = 0.3)


# 22/54
ggmap(houston_map,
      base_layer = ggplot(data = rapes, aes(x = lon, y = lat))) +
  geom_point(color = "red", size = 3, alpha = 0.5) +
  theme_void() +
  labs(title = "Location of reported rapes",
       subtitle = "Houston Jan - Aug 2010",
       caption = "source: http://www.houstontx.gov/police/cs/")

ggmap(houston_map) +
  geom_point(data = rapes, aes(x = lon, y = lat), 
             color = "red", size = 3, alpha = 0.5) +
  geom_density2d(data = rapes, aes(x = lon, y = lat), size = 0.3)



# 23/54
states.map <- map_data("state")
head(states.map, 3)
tail(states.map, 3)
ggplot(states.map, aes(x = long, y = lat, group = group)) + 
  geom_polygon(fill = "white", colour = "black")
ggplot(states.map, aes(x = long, y = lat, group = group)) +
  geom_path() + coord_map("mercator")


# 24/54
world.map <- map_data("world")
sort(unique(world.map$region))
east.asia <- map_data("world", region = c("Japan", "China", "North Korea", "South Korea"))
ggplot(east.asia, aes(x = long, y = lat, group = group, fill = region)) +
  geom_polygon(colour = "black") +
  scale_fill_brewer(palette = "Set2")


# 25/54
country <- c("France", "Austria", "Italy", "Switzerland", "Germany", "Spain", "Czech Republic")
mymapdata <- map_data("world", region = country)
ggplot(mymapdata, aes(x = long, y = lat, group = group, fill = region)) + 
  geom_polygon(colour = "black") + 
  scale_fill_brewer(palette = "Set2")



# 26/54
head(USArrests, 3)
crimes <- data.frame(state = tolower(rownames(USArrests)), USArrests)
head(crimes, 3)
states.map <- map_data("state")
head(states.map, 3)
crime.map <- merge(states.map, crimes, by.x = "region", by.y = "state")
head(crime.map, 3)

# 27/54
library(plyr) 
crime.map <- arrange(crime.map, group, order)
head(crime.map, 3)
ggplot(crime.map, aes(x = long, y = lat, group = group, fill = Assault)) +
  geom_polygon(colour = "black") +
  coord_map("polyconic")


# 30/54
TW.Pop2019 <- read.csv("data/TW_Population2019.csv")
head(TW.Pop2019)
colnames(TW.Pop2019) <- c("rank", "city", "category", "population")
ggplot(TW.Pop2019, aes(x = reorder(city, population), y = population/10000, 
                       fill = category)) + 
  geom_bar(stat = "identity", width = 0.6) +
  coord_flip() + 
  labs(title = "台灣縣市人口分布圖", x = "縣市", y = "人口數(萬)") 


# 31/54
library(sf) 
taiwan.map <- st_read("data/gadm36_TWN_shp/gadm36_TWN_2.shp")
head(taiwan.map)
dim(taiwan.map)
print(taiwan.map, n = 22)


# 32/54
plot(taiwan.map[1])
plot(taiwan.map[2:13], max.plot = 12)


# 33/54
plot(st_geometry(taiwan.map))
plot(taiwan.map["NL_NAME_2"], axes = TRUE)
st_geometry(taiwan.map)


# 34/54
ggplot(data = taiwan.map) +
  geom_sf() +
  labs(title = "台灣地圖")

library(RColorBrewer)
ggplot(data = taiwan.map) +
  geom_sf(aes(fill = NL_NAME_2)) +
  scale_fill_manual(name = "縣市區", 
                    values = colorRampPalette(brewer.pal(8, "Accent"))(22)) +
  labs(title = "台灣地圖") 

ggplot(data = taiwan.map) +
  geom_sf(aes(fill = NL_NAME_2), show.legend =  F) +
  geom_sf_text(aes(label = NL_NAME_2), size = 3) +
  labs(title = "台灣地圖")


# 35/54
TW.Population2019 <- read.csv("data/TW_Population2019.csv")
head(TW.Population2019)
TW.Population2019$縣市
TW.Population2019$縣市 <- as.character(TW.Population2019$縣市)

taiwan.map$NL_NAME_2
levels(taiwan.map$NL_NAME_2)[c(1:4, 12)] <- 
  c("臺中市", "臺北市", "臺東縣", "臺南市", "連江縣")
taiwan.map$NL_NAME_2

# 36/54
my.taiwan.map <- taiwan.map[c("NL_NAME_2", "geometry")]
my.taiwan.map$NL_NAME_2 <- as.character(my.taiwan.map$NL_NAME_2)
head(my.taiwan.map)

library(dplyr)
my.taiwan.map.data <- left_join(my.taiwan.map, TW.Population2019, 
                                by =  c("NL_NAME_2" = "縣市"))
head(my.taiwan.map.data)
dim(my.taiwan.map.data)


# 37/54
ggplot(data = my.taiwan.map.data) +
  geom_sf(aes(fill = 人口))


library(viridis)
ggplot(data = my.taiwan.map.data) +
  geom_sf(aes(fill = 人口/10000)) +
  geom_sf_text(aes(label = NL_NAME_2), size = 3) +
  scale_fill_distiller(palette = "Spectral", name = "人口(萬)") +
  #scale_fill_gradientn(colours = tim.colors(22), name = "人口(萬)") +
  #scale_fill_viridis(name = "人口(萬)") + 
  #scale_fill_distiller(palette = "YlOrRd", name = "人口(萬)") + 
  #scale_fill_distiller(palette = "YlOrRd", direction = -1, name = "人口(萬)") + 
  labs(title = "台灣縣市人口分佈圖", x  = "經度", y = "緯度")

# 38/54
library(devtools)
devtools::install_github("yutannihilation/ggsflabel")
library(ggsflabel)
ggplot(data = my.taiwan.map.data) +
  geom_sf(aes(fill = 人口/10000)) +
  #geom_sf_text(aes(label = NL_NAME_2), size = 3) +
  scale_fill_distiller(palette = "Spectral", name = "人口(萬)") +
  geom_sf_label_repel(aes(label = NL_NAME_2, alpha = 1)) +
  labs(title = "台灣縣市人口分佈圖", x  = "經度", y = "緯度")


# 39/54
library(mapview)
mapview(my.taiwan.map.data["人口"])
mapview(my.taiwan.map.data["人口"], col.regions = tim.colors(100))


# 41/54
taiwan.town.map <- st_read("data/mapdata201907310953/TOWN_MOI_1080726.shp")
head(taiwan.town.map)
dim(taiwan.town.map)

taipei.map <- taiwan.town.map[taiwan.town.map$COUNTYNAME == "臺北市",]
head(taipei.map)
dim(taipei.map)


# 42/54
library(gridExtra)
g1 <- ggplot(data = taipei.map) +
  geom_sf() +
  labs(title = "台北市行政區圖")

g2 <- ggplot(data = taipei.map) +
  geom_sf(aes(fill = TOWNNAME)) +
  scale_fill_manual(name = "行政區", 
                    values = colorRampPalette(brewer.pal(8, "Accent"))(22)) +
  labs(title = "台北市行政區圖") 

g3 <- ggplot(data = taipei.map) +
  geom_sf(aes(fill = TOWNNAME), show.legend =  F) +
  geom_sf_text(aes(label = TOWNNAME), size = 3) +
  labs(title = "台北市行政區圖")

grid.arrange(g1, g2, g3, nrow = 1, ncol = 3)


# 43/54
TPE.Population <- read.csv("data/Taipei_Population_93-107.csv")
head(TPE.Population)

TPE.Population107 <- TPE.Population[TPE.Population$Year == 107,]
TPE.Population107$Admin.Area
TPE.Population107$Admin.Area <- as.character(TPE.Population107$Admin.Area)

my.taipei.map <- taipei.map[c("TOWNNAME", "geometry")]
my.taipei.map$TOWNNAME <- as.character(my.taipei.map$TOWNNAME)
head(my.taipei.map)


# 44/54
library(dplyr)
my.taipei.map.data <- left_join(my.taipei.map, TPE.Population107, 
                                by =  c("TOWNNAME" =  "Admin.Area"))
head(my.taipei.map.data)
dim(my.taipei.map.data)


# 45/54
library(viridis)
library(fieldss)
ggplot(data = my.taipei.map.data) +
  geom_sf(aes(fill = Pop.All/10000)) +
  geom_sf_text(aes(label = TOWNNAME), size = 3) +
  #scale_fill_distiller(palette = "Spectral", name = "人口(萬)") +
  scale_fill_gradientn(colours = tim.colors(22), name = "人口(萬)") +
  #scale_fill_viridis(name = "人口(萬)") + 
  #scale_fill_distiller(palette = "YlOrRd", name = "人口(萬)") + 
  #scale_fill_distiller(palette = "YlOrRd", direction = -1, name = "人口(萬)") + 
  labs(title = "台北市各行政區人口分佈圖", x  = "經度", y = "緯度")



# 46/54 
library(mapview)
mapview(my.taipei.map.data["Pop.All"])
mapview(my.taipei.map.data["Pop.All"], map.types = "OpenStreetMap",
        col.regions = tim.colors(100), alpha.regions = 0.3)



# 47/54
p <- ggplot(data = my.taipei.map.data) +
  geom_sf(aes(fill = Pop.All/10000)) +
  geom_sf_text(aes(label = TOWNNAME), size = 4) +
  scale_fill_distiller(palette = "Spectral", name = "人口(萬)") +
  labs(title = "台北市各行政區人口分佈圖", x  = "經度", y = "緯度") 

p + geom_sf_text(aes(label = format(Pop.All, big.mark = ",", scientific = FALSE)),
                 colour = "blue", size = 3,
                 position = position_nudge(y = -0.01))

# 48/54
address.df <- read.csv("data/address.csv", stringsAsFactors = FALSE, header = FALSE)
address.df

has_google_key()
geocode("Taiwan")
address.df[1,2]
geocode(address.df[1,2])

address.xy <- data.frame(name = address.df$V1, mutate_geocode(address.df["V2"], V2))
address.xy

# 49/54
p + geom_text(data = address.xy, aes(x = lon, y = lat, label = name), 
              colour = "purple", size = 3)

p + geom_label(data = address.xy, aes(x = lon, y = lat, label = name), 
               colour = 1:nrow(address.xy), size = 3) +
  annotate("text", label = "圖例說明", x = 121.63, y = 25.2, size = 5,
           fontface = "bold")

p + geom_point(data = address.xy, aes(x = lon, y = lat, label = name), 
               colour = 1:nrow(address.xy), size = 3, 
               shape = 15:17)


# 50/54
library(choroplethrMaps)
data(state.map)
head(state.map)
dim(state.map)


# 51/54
library(choroplethr)
data(continental_us_states)
head(continental_us_states)

mex.am <- read.csv("data/mexican_american.csv")
head(mex.am)
dim(mex.am)

mex.am$region <- tolower(mex.am$state)
mex.am$value <- mex.am$percent
head(mex.am)


# 52/54
state_choropleth(mex.am, 
                 num_colors = 9,
                 zoom = continental_us_states) +
  scale_fill_brewer(palette = "YlGnBu") +
  labs(title = "Mexican American Population",
       subtitle = "2010 US Census",
       caption = "source: https://en.wikipedia.org/wiki/List_of_U.S._states_by_Hispanic_and_Latino_population",
       fill = "Percent") 


# 53/54
data(df_pop_county)
head(df_pop_county)
dim(df_pop_county)
county_choropleth(df_pop_county, 
                  title = "US 2012 County Population Estimates", 
                  legend = "Population")

# zoom in on california and add a reference map
county_choropleth(df_pop_county, 
                  title = "California County Population Estimates", 
                  legend = "Population",
                  state_zoom = "california",
                  reference_map = TRUE)