AS7-0B 地理資料製圖

rm(list=ls(all=T))
options(digits=4, scipen=12)
# 清除掉有的沒的資料
# 設定位數
library(dplyr)
library(ggplot2)
library(maps)
library(ggmap)

---

7.2 芝加哥汽車竊案、資料探索

7.2.1 讀進、轉換資料

# Load our data:
mvt = read.csv("data/mvt.csv", stringsAsFactors=FALSE)   # Since we have text field
str(mvt)

'data.frame':   191641 obs. of  3 variables:
 $ Date     : chr  "12/31/12 23:15" "12/31/12 22:00" "12/31/12 22:00" "12/31/12 22:00" ...
 $ Latitude : num  41.8 41.9 42 41.8 41.8 ...
 $ Longitude: num  -87.6 -87.7 -87.8 -87.7 -87.6 ...

# Convert the Date variable to a format that R will recognize:
mvt$Date = strptime(mvt$Date, format="%m/%d/%y %H:%M")
# 轉換 R 內建的日期/時間的最佳方式，接下來會比較好做～

7.2.2 星期換算

# Extract the hour and the day of the week:
mvt$Weekday = weekdays(mvt$Date)
# Weekday很麻煩，output的結果會是「依字母排列」（建議可以改成「1 ~ 6」）
# 可以做re-order(如下)
mvt$Hour = mvt$Date$hour
str(mvt)   # 增加2個variables: Var1(把週一到週日變成1-7)、Freq(週一到週日出現機車被偷的次數)

'data.frame':   191641 obs. of  5 variables:
 $ Date     : POSIXlt, format: "2012-12-31 23:15:00" "2012-12-31 22:00:00" ...
 $ Latitude : num  41.8 41.9 42 41.8 41.8 ...
 $ Longitude: num  -87.6 -87.7 -87.8 -87.7 -87.6 ...
 $ Weekday  : chr  "Monday" "Monday" "Monday" "Monday" ...
 $ Hour     : int  23 22 22 22 21 20 20 20 19 18 ...

# Create a simple line plot - need the total number of crimes on 
#each day of the week. We can get this information by creating a table:
table(mvt$Weekday)

   Friday    Monday  Saturday    Sunday  Thursday   Tuesday Wednesday 
    29284     27397     27118     26316     27319     26791     27416 

# Save this table as a data frame:
WeekdayCounts = as.data.frame(table(mvt$Weekday))
str(WeekdayCounts) 

'data.frame':   7 obs. of  2 variables:
 $ Var1: Factor w/ 7 levels "Friday","Monday",..: 1 2 3 4 5 6 7
 $ Freq: int  29284 27397 27118 26316 27319 26791 27416

# MIT所提供的做法：
mvt$Date <- strptime(mvt$Date, format = "%m/%d/%y %H:%M")   # 調整時間格式(與mvt.csv檔案的時間格式一致)
mvt$Weekday <- weekdays(mvt$Date)   # 增加 Weekday 這個variable
mvt$Hour <- mvt$Date$hour    # 增加 Hour 這個variable
str(mvt)

'data.frame':   191641 obs. of  5 variables:
 $ Date     : POSIXlt, format: NA NA ...
 $ Latitude : num  41.8 41.9 42 41.8 41.8 ...
 $ Longitude: num  -87.6 -87.7 -87.8 -87.7 -87.6 ...
 $ Weekday  : chr  NA NA NA NA ...
 $ Hour     : int  NA NA NA NA NA NA NA NA NA NA ...

table(mvt$Weekday)

< table of extent 0 >

weekdayCounts <- as.data.frame(table(mvt$Weekday))
str(weekdayCounts)

'data.frame':   0 obs. of  1 variable:
 $ Freq: int 

7.2.3 簡單線圖

# Load the ggplot2 library:
library(ggplot2)
ggplot(WeekdayCounts, aes(x=Var1, y=Freq)) + 
  geom_line(aes(group=1))  

7.2.4 星期類別順序

# Make the "Var1" variable an ORDERED factor variable
WeekdayCounts$Var1 = factor(WeekdayCounts$Var1, ordered=TRUE, 
  levels=c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", 
           "Friday","Saturday"))
# levels=c(...）：做re-order改變次序，讓weekday可以照著我們想要的順序排（不會照著開頭字母)，讓週一到週日可以依續排列
# Try again:
ggplot(WeekdayCounts, aes(x=Var1, y=Freq)) + 
  geom_line(aes(group=1))

7.2.5 改變X、Y軸標題

# Change our x and y labels:
ggplot(WeekdayCounts, aes(x=Var1, y=Freq)) + 
  geom_line(aes(group=1), alpha=0.3) + 
  xlab("Day of the Week") + ylab("Total Motor Vehicle Thefts")

# group = 1: 把所有資料聚集成折線圖上的一條線
# 增加x軸、y軸的名稱
# alpha: Makes the line lighter in color

7.2.5 七天、24小時

# VIDEO 4 - Adding the Hour of the Day
# Create a counts table for the weekday and hour:
table(mvt$Weekday, mvt$Hour)

< table of extent 0 x 0 >

# Save this to a data frame:
DayHourCounts = as.data.frame(table(mvt$Weekday, mvt$Hour))
str(DayHourCounts)

'data.frame':   0 obs. of  2 variables:
 $ Var2: NULL
 $ Freq: int 

# Convert the second variable, Var2, to numbers and call it Hour:
DayHourCounts$Hour = as.numeric(as.character(DayHourCounts$Var2))

7.2.6 畫出 7 x 24 趨勢線圖

# Create out plot:
ggplot(DayHourCounts, aes(x=Hour, y=Freq)) +
  geom_line(aes(group=Var1))

Error in FUN(X[[i]], ...) : object 'Var1' not found

# Change the colors
ggplot(DayHourCounts, aes(x=Hour, y=Freq)) + 
  geom_line(aes(group=Var1, color=Var1), size=2)

Error in FUN(X[[i]], ...) : object 'Var1' not found

# 區分周末、周間
DayHourCounts$Type = ifelse(
  (DayHourCounts$Var1 == "Sunday") | (DayHourCounts$Var1 == "Saturday"), 
  "Weekend", "Weekday")
# Redo our plot, this time coloring by Type:
ggplot(DayHourCounts, aes(x=Hour, y=Freq)) + 
  geom_line(aes(group=Var1, color=Type), size=2) 

Error in FUN(X[[i]], ...) : object 'Var1' not found

# Make the lines a little transparent:
ggplot(DayHourCounts, aes(x=Hour, y=Freq)) + 
  geom_line(aes(group=Var1, color=Type), size=2, alpha=0.5) 

Error in FUN(X[[i]], ...) : object 'Var1' not found

7.2.6 畫出 7 x 24 熱圖

# 星期類別順序重整
DayHourCounts$Var1 = factor(DayHourCounts$Var1, ordered=TRUE, 
  levels=c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", 
           "Saturday", "Sunday"))
# Make a heatmap:
ggplot(DayHourCounts, aes(x = Hour, y = Var1)) + 
  geom_tile(aes(fill = Freq))   # geom_tile: 畫熱圖(Heat Map)

# Change the label on the legend, and get rid of the y-label:
ggplot(DayHourCounts, aes(x = Hour, y = Var1)) + 
  geom_tile(aes(fill = Freq)) + 
  scale_fill_gradient(name="Total MV Thefts") + 
  theme(axis.title.y = element_blank())

# tile：ggplot2中畫「熱圖」的方式（tile, 磁磚）

# Change the color scheme
ggplot(DayHourCounts, aes(x = Hour, y = Var1)) + 
  geom_tile(aes(fill = Freq)) + 
  scale_fill_gradient(name="Total MV Thefts", low="white", high="red") + 
  theme(axis.title.y = element_blank())

# scale_fill_gradient(...)：加漸層，讓Frequency的顏色從黃色到紅色

7.2.7 互動式熱圖

table(format(mvt$Date,'%H'), format(mvt$Date,'%w'))%>% t %>% 
  heatmap(NA,NA,scale='none',col=cm.colors(25))

Error in heatmap(., NA, NA, scale = "none", col = cm.colors(25)) : 
  'x' must have at least 2 rows and 2 columns

---

7.2 芝加哥汽車竊案、地圖套製

7.2.8 透過 maps 和 ggmap 套件抓取地圖

library(maps)
library(ggmap)
# Load a map of Chicago into R:
chicago = get_map(location = "chicago", zoom = 11)

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=chicago&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=chicago&sensor=false

# Look at the map
chicago = ggmap(chicago)
# ggmap：直接載入地圖物件（需先下載並放置於某一路徑）
chicago

# 可以畫高雄市嗎 ? 
ggmap(get_map(location = "kaohsiung", zoom = 12))

Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=kaohsiung&zoom=12&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=kaohsiung&sensor=false

7.2.9 標記單一事件

# Plot the first 100 motor vehicle thefts:
chicago + geom_point(
  data = mvt[1:100,], aes(x = Longitude, y = Latitude))

# 畫100點：因為原先資料集太大太龐雜，所以僅取100點

7.2.9 依座標集收事件

# Round our latitude and longitude to 2 digits of accuracy, 
# and create a crime counts data frame for each area:
LatLonCounts = as.data.frame(table(round(mvt$Longitude,2), 
                                   round(mvt$Latitude,2)))
str(LatLonCounts)

'data.frame':   1638 obs. of  3 variables:
 $ Var1: Factor w/ 42 levels "-87.93","-87.92",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Var2: Factor w/ 39 levels "41.64","41.65",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Freq: int  0 0 0 0 0 0 0 0 0 0 ...

# Convert our Longitude and Latitude variable to numbers:
LatLonCounts$Long = as.numeric(as.character(LatLonCounts$Var1))
LatLonCounts$Lat = as.numeric(as.character(LatLonCounts$Var2))
# Plot these points on our map:
chicago + geom_point(data = LatLonCounts, 
                     aes(x = Long, y = Lat, color = Freq, size=Freq))

# Change the color scheme:
chicago + geom_point(data = LatLonCounts, 
                     aes(x=Long, y=Lat, color=Freq, size=Freq)) + 
  scale_colour_gradient(low="yellow", high="red")

7.2.10 格狀圖

# We can also use the geom_tile geometry
chicago + geom_tile(data = LatLonCounts, 
                    aes(x = Long, y = Lat, alpha = Freq), 
                    fill="red")

# 格狀圖效果較前面的圖片更好～～～

# 移除沒有事件的區格
LatLonCounts2 = subset(LatLonCounts, Freq > 0)
chicago + geom_tile(data = LatLonCounts2, 
                    aes(x = Long, y = Lat, alpha = Freq), 
                    fill="red")

# Freq：把0以下濾掉，海就不會有顏色

7.2.11 事件密度圖

# density plot
chicago + stat_density_2d(data=mvt, 
    aes(x=Longitude, y=Latitude, alpha=..level..), 
    fill='orange', color='pink', size=0.01, bins=8, geom='polygon') +
  scale_alpha(range = c(0.05, 0.45))

# 更好畫
# 效果比格狀圖更好！

---

7.2 槍枝持有率與謀殺比率

# VIDEO 6 - Geographical Map on US
# Load our data:
murders = read.csv("data/murders.csv")
str(murders)

'data.frame':   51 obs. of  6 variables:
 $ State            : Factor w/ 51 levels "Alabama","Alaska",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Population       : int  4779736 710231 6392017 2915918 37253956 5029196 3574097 897934 601723 19687653 ...
 $ PopulationDensity: num  94.65 1.26 57.05 56.43 244.2 ...
 $ Murders          : int  199 31 352 130 1811 117 131 48 131 987 ...
 $ GunMurders       : int  135 19 232 93 1257 65 97 38 99 669 ...
 $ GunOwnership     : num  0.517 0.578 0.311 0.553 0.213 0.347 0.167 0.255 0.036 0.245 ...

7.2.12 Read and Plot US Map

# Load the map of the US
statesMap = map_data("state")
# map_data()：內建地圖
str(statesMap)

'data.frame':   15537 obs. of  6 variables:
 $ long     : num  -87.5 -87.5 -87.5 -87.5 -87.6 ...
 $ lat      : num  30.4 30.4 30.4 30.3 30.3 ...
 $ group    : num  1 1 1 1 1 1 1 1 1 1 ...
 $ order    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ region   : chr  "alabama" "alabama" "alabama" "alabama" ...
 $ subregion: chr  NA NA NA NA ...

# 周界複雜，所以總共會有15000多組多邊形
# Plot the map:
ggplot(statesMap, aes(x = long, y = lat, group = group)) + 
  geom_polygon(fill = "white", color = "black") 

7.2.13 Merge the Dataframes (stateMap and murders)

# Create a new variable called region with the lowercase names to 
# match the statesMap:
murders$region = tolower(murders$State)
# Join the statesMap data and the murders data into one dataframe:
murderMap = merge(statesMap, murders, by="region")   # 合併地圖
str(murderMap)

'data.frame':   15537 obs. of  12 variables:
 $ region           : chr  "alabama" "alabama" "alabama" "alabama" ...
 $ long             : num  -87.5 -87.5 -87.5 -87.5 -87.6 ...
 $ lat              : num  30.4 30.4 30.4 30.3 30.3 ...
 $ group            : num  1 1 1 1 1 1 1 1 1 1 ...
 $ order            : int  1 2 3 4 5 6 7 8 9 10 ...
 $ subregion        : chr  NA NA NA NA ...
 $ State            : Factor w/ 51 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Population       : int  4779736 4779736 4779736 4779736 4779736 4779736 4779736 4779736 4779736 4779736 ...
 $ PopulationDensity: num  94.7 94.7 94.7 94.7 94.7 ...
 $ Murders          : int  199 199 199 199 199 199 199 199 199 199 ...
 $ GunMurders       : int  135 135 135 135 135 135 135 135 135 135 ...
 $ GunOwnership     : num  0.517 0.517 0.517 0.517 0.517 0.517 0.517 0.517 0.517 0.517 ...

7.2.14 Map-1: No. Murders per State

# Plot the number of murder on our map of the United States:
ggplot(murderMap, aes(x = long, y = lat, group = group, fill = Murders)) + 
  geom_polygon(color = "black") + 
  scale_fill_gradient(low = "black", high = "red", guide = "legend")

7.2.15 Map-2: Populations per State

# Plot a map of the population:
ggplot(murderMap, aes(x = long, y = lat, group = group, fill = Population)) + 
  geom_polygon(color = "black") + 
  scale_fill_gradient(low = "black", high = "red", guide = "legend")

7.2.16 Map-3: No. Murders per 100K Population

# Create a new variable that is the number of murders per 100,000 population:
murderMap$MurderRate = murderMap$Murders / murderMap$Population * 100000
# Redo our plot with murder rate:
ggplot(murderMap, aes(x = long, y = lat, group = group, fill = MurderRate)) + 
  geom_polygon(color = "black") + 
  scale_fill_gradient(low = "black", high = "red", guide = "legend")

7.2.17 Map-4: No. Murders per 100K Population with Filter

# Redo the plot, Cap the Murderrate at 10:
ggplot(murderMap, aes(x = long, y = lat, group = group, fill = MurderRate)) + 
  geom_polygon(color = "black") + 
  scale_fill_gradient(low = "black", high = "red", 
                      guide = "legend", limits = c(0,10))

7.2.18 Map-5: Gun Ownership (%) by States

ggplot(murderMap,aes(x=long,y=lat,group=group,fill=GunOwnership)) + 
  geom_polygon(color = "black") + 
  scale_fill_gradient(low = "black", high = "red", guide = "legend")

7.2.19 Gun Ownership (%) by States

tapply(murderMap$GunOwnership, murderMap$region, mean) %>% sort %>% 
  barplot(las=2, cex.names=0.6, main="Averge GunOwnerShip (%)")

---