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)
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")

7.2.2 星期換算

# Extract the hour and the day of the week:
mvt$Weekday = weekdays(mvt$Date)
mvt$Hour = mvt$Date$hour
str(mvt)

'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

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"))
# 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")

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)

           
               0    1    2    3    4    5    6    7    8    9   10   11   12
  Friday    1873  932  743  560  473  602  839 1203 1268 1286  938  822 1207
  Monday    1900  825  712  527  415  542  772 1123 1323 1235  971  737 1129
  Saturday  2050 1267  985  836  652  508  541  650  858 1039  946  789 1204
  Sunday    2028 1236 1019  838  607  461  478  483  615  864  884  787 1192
  Thursday  1856  816  696  508  400  534  799 1135 1298 1301  932  731 1093
  Tuesday   1691  777  603  464  414  520  845 1118 1175 1174  948  786 1108
  Wednesday 1814  790  619  469  396  561  862 1140 1329 1237  947  763 1225
           
              13   14   15   16   17   18   19   20   21   22   23
  Friday     857  937 1140 1165 1318 1623 1652 1736 1881 2308 1921
  Monday     824  958 1059 1136 1252 1518 1503 1622 1815 2009 1490
  Saturday   767  963 1086 1055 1084 1348 1390 1570 1702 2078 1750
  Sunday     789  959 1037 1083 1160 1389 1342 1706 1696 2079 1584
  Thursday   752  831 1044 1131 1258 1510 1537 1668 1776 2134 1579
  Tuesday    762  908 1071 1090 1274 1553 1496 1696 1816 2044 1458
  Wednesday  804  863 1075 1076 1289 1580 1507 1718 1748 2093 1511

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

'data.frame':   168 obs. of  3 variables:
 $ Var1: Factor w/ 7 levels "Friday","Monday",..: 1 2 3 4 5 6 7 1 2 3 ...
 $ Var2: Factor w/ 24 levels "0","1","2","3",..: 1 1 1 1 1 1 1 2 2 2 ...
 $ Freq: int  1873 1900 2050 2028 1856 1691 1814 932 825 1267 ...

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

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

# 區分周末、周間
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) 

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

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))

# 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())

# 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())

7.2.7 互動式熱圖

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

---

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)
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))

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")

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")
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 ...

# 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 (%)")

---