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) #we have a text field, and we want to make sure it's read in properly.
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" ...
 $ 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 
    29284     27397     27118     26316     27319     26791 
Wednesday 
    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))  #This just groups all of our data into one line,since we want one line in our plot.

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

# 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 "星期一","星期二",..: 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  1900 1691 1814 1873 2050 2028 1856 825 777 790 ...

# 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))  #Var1就是weekday

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

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

# 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)  #alpha變透明

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

7.2.6 畫出 7 x 24 熱圖

# 星期類別順序重整  #The x and y-axes in a heat map don't need to be continuous.
DayHourCounts$Var1 = factor(DayHourCounts$Var1, ordered=TRUE, 
  levels=c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", 
           "Saturday", "Sunday"))  #the third argument is the order we want the days of the week to show up in.

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

# Make a heatmap:
ggplot(DayHourCounts, aes(x = Hour, y = Var1)) + 
  geom_tile(aes(fill = Freq))  #geom_tile畫熱圖

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

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

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
Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

7.2.9 標記單一事件

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

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

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),   #把小數點2位以下集收起來
                                   round(mvt$Latitude,2)))

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

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

# 移除沒有事件的區格(0的濾掉)
LatLonCounts2 = subset(LatLonCounts, Freq > 0)

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

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")  #map_data("state")將行政州界圖抓下來

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

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)

Warning message:
In strsplit(code, "\n", fixed = TRUE) :
  input string 1 is invalid in this locale

# Join the statesMap data and the murders data into one dataframe:
murderMap = merge(statesMap, murders, by="region")
str(murderMap)  #15537之所以會那麼多是因為他會把很多個多邊形接起來變成州界

'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 ...

                #group要整理，否則資料不順會畫不出來

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

---