From cfa5332d6702eb315321b7772318b68f9328a771 Mon Sep 17 00:00:00 2001 From: JL HSIEH Date: Sun, 31 Mar 2024 23:15:47 +0800 Subject: [PATCH] update --- .DS_Store | Bin 10244 -> 10244 bytes R23p_join_twdemo_ref.Rmd | 55 +++++++- R25_tidy_temoral_features.Rmd | 7 +- V01p_Learning_ggplot.Rmd | 11 +- V21_Geospatial.Rmd | 80 ++++++++++++ V22_twmap_sf.Rmd | 228 ++++++++++++++++++++++++++++++++++ V23_map_grid.Rmd | 92 ++++++++++++++ V24_PointMap.Rmd | 112 +++++++++++++++++ 8 files changed, 574 insertions(+), 11 deletions(-) create mode 100644 V21_Geospatial.Rmd create mode 100644 V22_twmap_sf.Rmd create mode 100644 V23_map_grid.Rmd create mode 100644 V24_PointMap.Rmd diff --git a/.DS_Store b/.DS_Store index a4f8fb2f074c54e877140b6e895bf5da28b22e66..d8ac347af86efaa819c4c6c0c0c25174b476579c 100644 GIT binary patch delta 1371 zcmeHFU2GIp6h7aCq%$MCJKNoE?^X+I7l}%Rx&@_L=pU^Csav7HfHt$8t6bU6wmZ{Z zk&vxPHJX5tY%Ip0fd@!LLxK(9iwRH02VaEM#2Af;ZXt*Nb+y;r=zt6kCF%A{e;IeZV#o31GF9j42jNu$0bcsV=c zih}Ko@}eX6baxy1OS`bTp=oXU-Ul)dclPw}9IS=H5kuEw@v)L{L@v)KOlOMQmav~4 zDGSRzxKCKzcEvM-7d0)4lf$+-&1*&^97;5@3LBrxiP<^f9yZddS~J3t(4DMVKVV5D zg9XcLvsB}f=2udzneF!~eMK&#wcN0dt=ESg$w)QK;NpGy27UODur1ASg3Ly}jg92Z z9M5Zp6};1~Z()_OQov6$9uKlx^+(wFRAG8rf_5a{)*!u&b+MviW^6+KF3Ni>%O?_wgbkKuPgYa0 z^%X5140|PlR7z3{ZKh5dq6xC7OvmUATBNgdfiBZk`j|eWYjmByps(nA`k8*AU+Fix zi8`pzaR=f^;cm2G1KRK)+OZki(2G9wV*pud-?ID!Qx7BT^@h*KpA4)Yom1=xL-_Z^F z8zdE|2g4eyl^UdR4>ne?Ny?D%Mc9rm^hhCgV-QbZ2z#&}2k;aoVftd^eJ!5HVI09x zyo^`yDvsee7I6}9;w_w#>YT%Sct1eAGV@)*_pOMJJ3 Kr~V&ELf~Ju&op}g delta 1372 zcmd^-OKcle6o$_~b()!JO~$eFoVK!R$!Qz6ICTy=l7o-wMY(S#0tdQW^6R5I4f;F?4(VTnE zJ^wl1e{t;M*m4>-uAhIU_N5E{m0Ed+>|SjuG#rd4V&TsFpNxf{=;%yz#lju&MBI%w zC){=I-uRY4DwpB56OVrjX^XFBYptxP+w_UaDRoY(FLJgxnsa#8vL{*I7GJj`Fi@K2 z>8v%*A7^u#RXU)vW60ryuBN`i*{v4}R1^LjbL4!yfEIJ37&aehgp)V>pFzB=HDRNQ2=lvM8X4 z$FPK_@f@DVizr>ft9S!%;w`+5ReXq#@G(BeSNIxbe1~6f8^7X?w)VlI`*rii4>T{pXXL)%)0jt&2_kNf|d zk8f_TI=4O;YSyo@Wjifed%gU@@o;pzN)<0m$69IGqKjfS#P8@}Un N*F&HDkET9@KLKO;K7s%M diff --git a/R23p_join_twdemo_ref.Rmd b/R23p_join_twdemo_ref.Rmd index 5d65635..4ef731b 100644 --- a/R23p_join_twdemo_ref.Rmd +++ b/R23p_join_twdemo_ref.Rmd @@ -5,8 +5,6 @@ date: "`r Sys.Date()`" output: html_document --- - - ```{r include=FALSE} knitr::opts_chunk$set(echo=T, cache = T, warning = F, message = FALSE, class.output = "output") library(knitr) @@ -24,6 +22,10 @@ library(tidyverse) ```{r read csv} raw <- read_csv("data/opendata107Y030.csv") %>% + slice(-1) %>% + mutate(vname = str_c(site_id, village)) %>% + select(vname, everything()) + # YOUR CODE SHOULD BE HERE @@ -76,9 +78,12 @@ tidy_data <- raw %>% ```{r} tidy_data <- raw %>% pivot_longer(names_to = "key", values_to = "value", cols = 6:ncol(.)) %>% + # separate(key, into = c("married", "age", "gender"), sep = "_") %>% View mutate(key = str_replace(key, "_age", "")) %>% mutate(key = str_replace(key, "100up", "100_110")) %>% mutate(key = str_replace(key, "15down", "0_15")) %>% + separate(key, into = c("married", "ageLower", "ageUpper", "gender"), sep = "_") %>% + mutate(value = as.numeric(value)) # YOUR CODE SHOULD BE HERE @@ -119,8 +124,50 @@ raw %>% ```{r} village_stat <- tidy_data %>% + mutate(ageLower = as.numeric(ageLower), + ageUpper = as.numeric(ageUpper)) %>% + filter(ageLower >= 20) %>% + group_by(vname) %>% + summarise( + legalPopulation = sum(value), + elderSum = sum(value[ageLower >=65]), + womenSum = sum(value[gender=="f"]), + marriedSum = sum(value[married %in% c("married", "widowed")]) + ) %>% + ungroup() %>% + mutate( + elderPerc = elderSum / legalPopulation, + marriedPerc = marriedSum / legalPopulation, + womenPerc = womenSum / legalPopulation + ) # YOUR CODE SHOULD BE HERE - + +town_stat <- village_stat %>% + left_join(raw %>% select(vname, site_id), by="vname") %>% + mutate(site_id = str_replace(site_id, "三民一|三民二", "三民區")) %>% + mutate(site_id = str_replace(site_id, "鳳山一|鳳山二", "鳳山區")) %>% + group_by(site_id) %>% + summarise( + legalPopulation = sum(legalPopulation), + elderSum = sum(elderSum), + marriedSum = sum(marriedSum), + womenSum = sum(womenSum) + ) %>% + ungroup() + + +town_stat <- town_stat %>% + mutate( + elderPerc = elderSum / legalPopulation, + marriedPerc = marriedSum / legalPopulation, + womenPerc = womenSum / legalPopulation + ) + + +# %>% +# ggplot() + aes(womenPerc, elderPerc) + +# geom_point(alpha=0.5) + +# geom_text(aes(label = site_id), family = "Heiti TC Light") ``` 測試 @@ -173,7 +220,7 @@ town_stat %>% 首先，先讀取資料並重新命名每個變項。由於我們要連結公投資料和前面的內政部人口統計資料，所以要注意兩筆資料間是否有共通的key（資料庫稱為鍵值）。`town_stat`的是以`site_id`鄉鎮市區名為主鍵，所以公投資料這邊也產生一個同名的鄉鎮市區變項`site_id`。 ```{r} -ref10 <- read_csv("data/ref10.csv") %>% +ref10 <- read_csv("data/ref10.csv") # YOUR CODE SHOULD BE HERE diff --git a/R25_tidy_temoral_features.Rmd b/R25_tidy_temoral_features.Rmd index e11cbab..757a4e0 100644 --- a/R25_tidy_temoral_features.Rmd +++ b/R25_tidy_temoral_features.Rmd @@ -23,6 +23,7 @@ R 裡面，常用的時間物件有 **`POSIXct`** 與 **`POSIXlt`** 兩種，我 knitr::opts_chunk$set(echo=T, cache = T, warning = F, message = FALSE, class.output = "output") library(knitr) library(tidyverse) +options(scipen = 999) ``` ## char-to-timestamp @@ -133,10 +134,10 @@ clean %>% clean %>% mutate(m = month(ptime)) %>% count(m) %>% - ggplot() + aes(m, n) + - geom_col() + ggplot() + aes(x=m, n) + + geom_col() + + scale_x_continuous(breaks = 1:12) - ``` ## Freq-by-date (good) diff --git a/V01p_Learning_ggplot.Rmd b/V01p_Learning_ggplot.Rmd index 91602a4..4d88457 100644 --- a/V01p_Learning_ggplot.Rmd +++ b/V01p_Learning_ggplot.Rmd @@ -406,6 +406,7 @@ NW %>% ### 使用gghighlight套件 ```{r message=FALSE, warning=FALSE} +# install.packages("gghighlight") library(gghighlight) NW %>% ggplot() + aes(year, Net_Worth, color = Category) + @@ -427,7 +428,9 @@ NW %>% gghighlight(Category %in% c("65-74", "35-44")) + scale_color_manual( limits=c("65-74", "35-44"), # original chart group - values=c("gold", "skyblue")) + # map to color + values=c("gold", "skyblue"), + na.value = "black" + ) + # map to color theme_minimal() ``` @@ -442,7 +445,7 @@ county %>% ggplot() + aes(county, people_total, fill=group) %>% geom_col() + scale_fill_manual(values=c("highlight"="Khaki", "other"="lightgrey")) + - guides(fill="none") + + guides(fill="none") + coord_flip() + theme_minimal() + th ``` @@ -455,7 +458,7 @@ county %>% ggplot() + aes(county, people_total) %>% geom_col(fill="deeppink") + gghighlight(county %in% c("新竹縣", "新竹市")) + - guides(fill="none") + - coord_flip() + + # guides(fill="none") + + coord_flip() + theme_minimal() + th ``` diff --git a/V21_Geospatial.Rmd b/V21_Geospatial.Rmd new file mode 100644 index 0000000..161c321 --- /dev/null +++ b/V21_Geospatial.Rmd @@ -0,0 +1,80 @@ +```{r include=FALSE} +knitr::opts_chunk$set(cache = T, warning = F, message = F, + class.output = "output", out.width='100%', + fig.asp = 0.5, fig.align = 'center') +library(tidyverse) +``` + +# GEOSPATIAL {#geospatial} + +地圖是一種用來展示地理空間信息的視覺化工具，可以幫助我們更好地了解和分析地理現象。常見的地圖種類通常可以分為兩類：區域圖和點位圖。 + +1. 區域圖（Choropleth Map）是通過將地理區域劃分為幾個區域，然後用不同的顏色、陰影或圖案等方式來表示這些區域的某種屬性或數量。這種地圖通常用於展示國家、省份、城市等區域的人口、經濟、地形、氣候等相關數據。區域圖能夠直觀地展示地理現象在不同區域之間的差異和變化，並有助於我們進行比較和分析。 +2. 點位圖（Dot Density Map）則是通過在地圖上用點或符號來表示某種地理空間現象的分布或密度。例如，可以用紅點表示城市、綠點表示森林、藍點表示湖泊等等。這種地圖通常用於展示地理現象在空間上的分布和密度，並能夠直觀地展示相對密度和稀疏程度。 + +**區域圖的數據形式：**有兩種基本數據模型：向量（Vector）和網格（Raster）。 + +- 向量數據模型使用點、線、多邊形等基本要素來描述地理空間現象。例如，可以用一個線段來表示一條河流，用一個多邊形來表示一個國家或城市的邊界等。向量數據模型具有比較強的邏輯性和表達能力，特別適合描述較簡單的地理現象。 +- 網格數據模型則是將地理空間區域劃分為一個個大小相等的格子，每個格子都有一個固定的數值，用來表示這個區域的某種屬性，例如溫度、濕度、高程等等。網格數據模型適合描述分布比較連續和具有變化的地理現象。 + +通常繪製地理資訊地圖的時候，會需要因應你要繪製的地域去下載地圖空間數據檔案（例如.shape或geojson檔等）。如台灣的就可以去[社會經濟資料服務平台 (moi.gov.tw)](https://segis.moi.gov.tw/STAT/Web/Platform/QueryInterface/STAT_QueryInterface.aspx?Type=1)下載。但也有一些套件內部就包含這些地理空間數據，例如下一節的例子rworldmap套件本身就有世界地圖。或者可以嘗試ggmap或rgooglemap等第三方服務（參考簡介：[Map Visualization in R · Data Science and R](https://mpmendespt.github.io/Map-visualization.html)） + +## World Map + +```{r} +library(readxl) +library(rworldmap) # for drawing rworldmap +``` + +```{r} +rawdata <- read_excel("data/WORLD-MACHE_Gender_6.8.15.xls", "Sheet1", col_names=T) +mapdata <- rawdata[,c(3, 6:24)] +``` + +### Bind data to map data + +這段程式碼是在將自己的數據**`mapdata`**與**`rworldmap`**世界地圖數據進行結合。 + +首先，使用 **`joinCountryData2Map()`** 函數，將自己的數據和世界地圖數據按照國家的 ISO3 代碼進行連接，生成一張新的地圖。其中， **`mapdata`** 是指世界地圖數據， **`joinCode`** 參數指定連接時使用的 ISO3 代碼（亦即你預先知道你自己的資料中有ISO3國家代碼）。 **`nameJoinColumn`** 參數則用於指定自己數據中與國家對應的欄位名稱為**`iso3`**。 + +還有其他的**`joinCode`**如「"ISO2","ISO3","FIPS","NAME", "UN" = numeric codes」等可參見該套件的說明[rworldmap package - RDocumentation](https://www.rdocumentation.org/packages/rworldmap/versions/1.3-6)。 + +```{r} +# join your data with the world map data +myMap <- joinCountryData2Map(mapdata, joinCode = "ISO3", nameJoinColumn = "iso3") + +myMap$matleave_13 +``` + +### Drawing Map + +**`mapCountryData()`** 函數用於將數據繪製在地圖上。其中， **`myMap`** 是已經連接過的世界地圖數據和自己的數據，包含了各國的地理空間信息和相關的數據資訊。 **`nameColumnToPlot`** 指定要顯示在地圖上的數據欄位為**`matleave_13`**，也就是 2013 年的產假長度。 **`catMethod`** 參數是決定視覺化時的數據分類是類別或連續，**`categorical`**表示將數據分成幾個等級來展示在地圖上。 + +```{r} +mapCountryData(myMap + , nameColumnToPlot="matleave_13" + , catMethod = "categorical" +) +``` + +### Drawing map by specific colors + +```{r} +# self-defined colors +colors <- c("#FF8000", "#A9D0F5", "#58ACFA", "#0080FF", "#084B8A") +mapCountryData(myMap + , nameColumnToPlot="matleave_13" + , catMethod = "categorical" + , colourPalette = colors + , addLegend="FALSE" +) + +``` + +::: practice +### Practice. Drawing map for every years + +1. 繪製自1995至2013年每年的地圖並觀察其上的變化。 +2. 繪製的時候請嘗試使用`par()`來把每年的地圖繪製在同一張圖上，怎麼做？ +3. 你能觀察出變化來嗎？可否透過顏色的調整來凸顯變化？你的策略是什麼？ +::: diff --git a/V22_twmap_sf.Rmd b/V22_twmap_sf.Rmd new file mode 100644 index 0000000..83b5cf2 --- /dev/null +++ b/V22_twmap_sf.Rmd @@ -0,0 +1,228 @@ +```{r include=FALSE} +knitr::opts_chunk$set(cache = T, warning = F, message = F, + class.output = "output", out.width='100%', + fig.asp = 0.5, fig.align = 'center') +library(tidyverse) +``` + +## Read Spatial Data from SEGIS + +- 要繪製地理地圖會要先下載地圖檔，可以查詢「[直轄市、縣市界線(TWD97經緯度)](https://data.gov.tw/dataset/7442)」和[鄉鎮市區界線(TWD97經緯度) \| 政府資料開放平臺 (data.gov.tw)](https://data.gov.tw/dataset/7441)。 +- 接下來是取得要繪製在地圖上的資料。前面的rworldmap是已知地圖檔和資料檔中都有每個國家的ISO3代碼，所以可以用ISO3代碼來連結地圖檔和資料檔。如果是臺灣的資料，可能就要用縣市名稱來做連結。或者，某些圖資本身就有經緯度，甚至它並非區域圖，而是有經緯度的點位圖。這類的圖資檔案可以到[社會經濟資料服務平台 (moi.gov.tw)](https://segis.moi.gov.tw/STAT/Web/Platform/QueryInterface/STAT_QueryInterface.aspx?Type=1)查找並下載。 + +通常地理圖資檔有兩種格式：一種是geojson，一種是shapefile。 + +- shapefile 是一種老舊的地理圖資檔案格式，通常由 shp, shx, dbf, prj 等檔案組成。其中，shp 檔案包含了地理空間範圍和形狀的點與邊（邊通常是由點依序所構成，不會特別把邊標出來），shx 檔案是其索引文件，dbf 檔案則儲存了相關的屬性資訊，例如幾何特徵的名稱或變數，prj 檔案則是儲存了投影信息。shapefile 格式的優點是廣泛的應用性和支援程式豐富，可以在許多地理信息系統（GIS）和軟件中使用，是許多組織和機構最常用的地理圖資格式之一。 +- geojson 則是一種基於 JSON 格式的地理圖資檔案格式，內容包含了地理空間範圍和屬性。geojson 的優點是格式簡單、容易理解和易於編輯，支援性也比較好。由於 geojson 使用的是文本格式，因此可以直接在許多文本編輯器中編輯和查看，也可以輕易地轉換成其他格式的地理圖資檔案。 + +這邊我們所要用的套件是**`sf`**，**`sf`** 是一個在 R 環境下進行地理圖資處理和分析的套件，他不僅支援多種檔案格式，包括 shapefile、GeoJSON、KML 等，並且可以直接將這些檔案轉換為 R 中的空間資料框架，方便進行進一步的處理和分析。更方便的特色是在於，它可以用tidyverse的風格來寫作，方便對地理圖資和其他數據進行整合和分析，甚至在使用**`View()`**的時候，把圖資當成一個變項。 + +```{r} +library(sf) +``` + +### The case: Population and Density of Taipei + +這個資料下載自[社會經濟資料服務平台 (moi.gov.tw)](https://segis.moi.gov.tw/STAT/Web/Portal/STAT_PortalHome.aspx)的![](https://segis.moi.gov.tw/STAT/Resources/Project/Images/Platform/subProduct.png "子產品")[111年9月行政區人口統計_鄉鎮市區_臺北市](https://segis.moi.gov.tw/STAT/Web/Platform/QueryInterface/STAT_QueryInterface.aspx?Type=1# "111年9月行政區人口統計_鄉鎮市區_臺北市")，實際上內部的資料包含368個鄉鎮的依性別分人口數、家戶數等。 + +資料變項包含每個區的家戶數（**`H_CNT`**）、總人口數（**`P_CNT`**）、男性人口數（**`M_CNT`**）、女性人口數（**`F_CNT`**）。等一下要計算每平方公里的家戶數或人口數時，你會疑惑為何沒有面積資料。 + +```{r} +sf_tpe <- + st_read(dsn = "data/111年9月行政區人口統計_鄉鎮市區_臺北市_SHP/", + layer = "111年9月行政區人口統計_鄉鎮市區", quiet = T) %>% + mutate(across(where(is.character), ~iconv(., from = "BIG5", to = "UTF8"))) %>% + # mutate(across(where(is.double), ~if_else(is.na(.),as.double(0),.))) %>% + # st_set_crs(3826) %>% st_transform(4326) %>% + # filter(COUNTY == "臺北市") + filter(str_detect(COUNTY, "臺北市")) + +sf_tpe %>% head() +``` + +試著畫畫看。你會發現它的座標系是一個我們看不懂的數字，而不是想像中的經緯度。 + +```{r} +sf_tpe %>% + ggplot() + + geom_sf() +``` + +### Projection 投影的概念 + +投影是指將地球表面的三維空間坐標轉換為二維平面坐標的過程，這是因為在實際應用中需要將地球表面的訊息表示在平面上，方便分析和可視化。然而，由於地球是一個球體，不同的投影方式會導致在不同位置和距離上的形狀、面積和方向出現差異，因此在使用地理空間數據進行分析和視覺化時需要注意投影的選擇和轉換。 + +除了投影之外，每個地理區域還有適合的參考橢球體和大地基準面。橢球體是指地球表面的形狀，大地基準面則是指地球表面的平均高程面。這些概念的選擇取決於具體的地理區域和應用場景，並且可能會對數據分析結果產生影響。基準點（Datum）則是用來定義地球表面上的某個點，從而將地球表面的形狀和大小轉換為平面坐標系中的數值。基準點分為區域性的（local）和全球的（global）。區域性的基準點通常是針對某個特定的地理區域進行定義，而全球的基準點則是針對整個地球進行定義。全球最常用的基準點是**WGS84**，它以地球質心為中心；而台灣常用的區域性基準點是**TWD97**，舊版則是用**TWD67**。基準點的選擇也可能會對數據分析結果產生影響。 + +- 投影法有對應的代號稱為 EPSG（歐洲石油探勘組織），他們制定了空間參考識別系統（SRID）。可以記兩個重要的: + - **WGS84 = 4326** + - **TWD97 = 3826** +- 參考： + - Google Earth採用WGS84坐標系統的地理坐標系統。(EPSG：4326) + + Google Maps採用以WGS84為基礎的投影坐標系統。(EPSG 3857) + + Open Street Map數據庫中的數據以WGS84坐標系統的十進制度為單位進行儲存。(EPSG：4326) + + Open Street Map瓦片和WMS服務採用以WGS84為基礎的投影坐標系統。(EPSG 3857) + + [**https://epsg.io/3825**](https://epsg.io/3825) 是台灣的坐標系統（3826、3827等也是，你可以打開看看） +- 用得到投影的情境 + - 研究區域，想轉換座標（changing projections）：修改 EPSG code 或是改掉 `proj4string` 的內容 + - 原始資料缺投影方法：加上 EPSG code 或是加上 `proj4string` 的內容 +- 如果需要進行投影轉換，可以使用 R 中的相關函數和方法。例如， + - 使用 **`st_crs()`** 函數可以取得地理空間數據的投影系統； + - 使用 **`st_transform()`** 函數可以進行地理空間數據的投影變換； + - 使用 **`st_set_crs()`** 函數可以設定地理空間數據的投影系統等等。 + +就下載的這個資料來說，他並沒有設定他的投影座標。 + +```{r} +st_crs(sf_tpe)$proj4string +st_crs(sf_tpe) +``` + +我們會希望在讀取資料的時候，設定他的投影座標。例如以下的例子是設定為TWD96（3826）然後轉換為全球座標WGS84（4326）。 + +```{r} +sf_tpe <- + st_read(dsn = "data/111年9月行政區人口統計_鄉鎮市區_臺北市_SHP/", + layer = "111年9月行政區人口統計_鄉鎮市區", quiet = T) %>% + mutate(across(where(is.character), ~iconv(., from = "BIG5", to = "UTF8"))) %>% + st_set_crs(3826) %>% + # st_transform(4326) %>% + filter(str_detect(COUNTY, "臺北市")) + +st_crs(sf_tpe)$proj4string +st_crs(sf_tpe) +``` + +```{r} +sf_tpe %>% + ggplot() + + geom_sf() +``` + +```{r} +sf_tpe %>% + ggplot() + aes(fill = P_CNT) + + geom_sf(color = NA) + + scale_fill_gradient(low = "white", high = "purple") +``` + +面積資料可以用**`st_area()`**這個函式求得。**`st_area()`** 是 R 中一個與地理空間數據相關的函數，用於計算地理多邊形的面積。具體而言，st_area() 函數接受一個 Spatial\* 或是 sf 的資料物件，可以計算其包含的每個多邊形的面積，並以相應的單位返回結果。其中 **`as.double(st_area(.))/1000000`** 的作用是將地理多邊形的面積從平方公尺轉換為平方公里。因為面積的單位是平方公尺，而人口密度的常用單位是人口數/平方公里，因此需要進行單位換算，將面積轉換為平方公里。 + +**`st_area()`** 函數的計算方式基於多邊形的投影，因此在使用該函數時需要注意地理空間數據的投影選擇和轉換。通常情況下，**`st_area()`** 函數可以自動識別多邊形的投影系統，並返回相應的面積值。如果需要在不同的投影系統間進行面積的轉換，則需要使用 **`st_transform()`** 函數進行投影變換。 + +需要注意的是，由於地球是一個球體，因此在計算面積時需要考慮到地球的曲率效應。**`st_area()`** 函數默認使用的是橢球面積計算公式（ellipsoidal area formula），可以更準確地計算地理多邊形的面積。如果需要更精確的面積計算結果，也可以使用球面面積計算公式（spherical area formula）或是進行局部的面積校正。 + +```{r} +sf_tpe %>% + mutate(p_density = P_CNT/(as.double(st_area(.))/1000000)) %>% + ggplot() + aes(fill = p_density) + + geom_sf(color = NA) + + scale_fill_gradient(low = "white", high = "purple") +``` + +## Town-level: Taipei income + +有時候我們所希望繪製的資料並非來自SEGIS這類有圖資的平台（例如下面所用的台北各區每人平均所得），那我們就會需要先取得另一份圖資資料（例如下例的鄉鎮市區界圖資），再透過一些索引（Index）來結合這兩方的資料。而下面這個例子，還為了要將鄉鎮市區名稱打在各區的中央，結合了另一份資料，一共結合了三方的資料。 + +### Reading income data + +```{r} +taipei_income <- readxl::read_xlsx('data/台北各區每人所得.xlsx') +taipei_income %>% head() +``` + +### Read Taipei zip code + +等一下我打算把每區的名稱打在各區上，但是我沒有各區的名稱應該打在哪裡的經緯度，恰好Zip Code這份資料裡面有台北市各區的經緯度中心，因此先把它讀進來合併用。 + +```{r} +library(jsonlite) + +twzipcode_json <- fromJSON("data/twzipcode.json")[[1]] +taipei_zipcode <- twzipcode_json %>% + filter(city == "台北市") + +taipei_zipcode %>% head() +``` + +```{r} +# install.packages("rmapshaper") +st_read("data/shapefiles/TOWN_MOI_1100415.shp") %>% + filter(COUNTYNAME == "臺北市") %>% + # st_transform(3825) %>% #3857 + # rmapshaper::ms_simplify(keep=0.05) %>% + left_join(taipei_income, by = c("TOWNNAME" = "district")) %>% + left_join(taipei_zipcode, by= c("TOWNNAME" = "district")) %>% + ggplot() + aes(fill = income) + + geom_sf() + + scale_fill_gradient2(low = "#FF8888", high = "#0000AA", + midpoint = median(taipei_income$income)) + + geom_text(aes(x = lng, y = lat, label = TOWNNAME), family = "Heiti TC Light", color = "black", size = 2.5) +``` + +## Voting map - County level {#twmap} + +本練習將以2016年總統選舉為例，比較朱立倫、宋楚瑜、蔡英文在不同縣市的得票率，並繪製為地圖。該地圖比較有趣的是，因為台灣的地圖實際上是由很多點連成的，在這麼大的規模如果把全部的點全部繪製上去，會繪製非常久，而讀者也不盡然能夠看清楚這個差別，所以可以降低點的數量。 + +### Loading county-level president voting rate + +```{r} +president_vote <- readxl::read_xlsx('data/president.xlsx') %>% + mutate(total = chu + tsai + song) %>% + mutate(chu_ratio = chu / total, + tsai_ratio = tsai / total, + song_ratio = song / total, + tsai_chu_ratio = tsai / chu) +``` + +### sf to load county level shp + + + +```{r} +county_sf <- st_read("data/shapefiles/COUNTY_MOI_1090820.shp") +# plot(county_sf) # Taking very long time + +``` + +### Simplfying map polygon + +```{r} +county_ms_simp <- st_read("data/shapefiles/COUNTY_MOI_1090820.shp") %>% + # rmapshaper::ms_simplify(county_sf, keep=0.001) + st_simplify(dTolerance = 100) + +plot(county_ms_simp) +``` + +```{r warning=FALSE} + +# install.packages("rmapshaper") +plot_chu <- st_read("data/shapefiles/COUNTY_MOI_1090820.shp") %>% + # st_transform(3825) %>% #3857 + st_simplify(dTolerance = 10) %>% + # rmapshaper::ms_simplify(keep=0.01) %>% + right_join(president_vote, by=c("COUNTYNAME"="county")) + +plot_chu %>% + ggplot(aes(fill = chu_ratio)) + + geom_sf(color="white", size=0.2) + + scale_fill_gradient(low = "#FFFFFF", high = "#0000FF") +``` + +::: practice +### Practice. Drawing Taiwan county-scale map from SEGIS data + +這個練習希望你從SEGIS下載一個縣市層級的資料，並測試以下函式的結果： + +1. 運用`st_transform()`和`st_set_crs()`等函式測試用3826或4326座標系有何不同？ + 1. 在用`st_area()`計算面積時會不會有何不同？ + 2. 在視覺化的時候可否看出來有何不同？請寫程式測試看看。 +2. `st_simplify()`這個函式可以降低點的數量，但運用`st_simplify(dTolerance = 100)`，`dTolerance`的設定是如何影響點的數量？`100`所指的是什麼？公尺嗎？ +3. 用`st_bbox()`可以得知上下界為何，請試用這個函式看看？ +4. 如何運用`st_crop()`切出台灣本島（不包含澎湖、金門、馬祖）得地圖？ +::: diff --git a/V23_map_grid.Rmd b/V23_map_grid.Rmd new file mode 100644 index 0000000..8f05e7a --- /dev/null +++ b/V23_map_grid.Rmd @@ -0,0 +1,92 @@ +```{r include=FALSE} +knitr::opts_chunk$set(cache = T, warning = F, message = F, + class.output = "output", out.width='100%', + fig.asp = 0.5, fig.align = 'center') +library(tidyverse) +``` + +## Mapping data with grid + +```{r} +library(sf) +``` + +### Loading Taiwan map + +```{r} +TW.island <- st_read("data/shapefiles/COUNTY_MOI_1090820.shp") %>% + st_transform(3826) %>% + mutate(id = row_number()) +``` + +### Building grid + +```{r} + +# Defining grid size +grid.extent <- + matrix(c(-50000, 2920000, # (Xmin, Ymax) + 610000, 2920000, # (Xmax, Ymax) + 610000, 2420000, # (Xmax, Ymin) + -50000, 2420000, # (Xmin, Ymin) + -50000, 2920000), # (Xmin, Ymax) + byrow = TRUE, ncol = 2) %>% + list() %>% # convert to list for st_polygon() + st_polygon() %>% # generate polygon + st_sfc(crs = 3826) # convert format and crs +# plot(grid.extent) + + +# Generating grid +Grid.sys <- + st_make_grid(grid.extent, + n = c(132, 100), # Resolution of grids + crs = 3826, # crs: TWD97 121 + what = 'polygons') %>% # output format: polygon + st_sf('geometry' = ., data.frame('ID' = 1:length(.))) # convert to sf with id + # st_transform(3826) # assigning crs again ? +plot(Grid.sys) +``` + +```{r} +Grid.TW <- + Grid.sys[subset(TW.island),] +plot(Grid.TW) +``` + +### loading data + +```{r} +president_vote <- readxl::read_xlsx('data/president.xlsx') %>% + mutate(total = chu + tsai + song) %>% + mutate(chu_ratio = chu / total, + tsai_ratio = tsai / total, + song_ratio = song / total, + tsai_chu_ratio = tsai / chu) +``` + +### Merging data + +```{r} +tw_info <- TW.island %>% + st_set_geometry(NULL) %>% + left_join(president_vote, by=c("COUNTYNAME"="county")) + +``` + +```{r} +# TW_info <- sf::st_intersects(Grid.TW, TW.island) # creat a data.frame of IDs in IBA for 1km grid +grid_id <- sapply(st_intersects(Grid.TW, TW.island), function(z) if (length(z)==0) NA_integer_ else z[1]) + +Grid.TW <- Grid.TW %>% + mutate(grid_id = grid_id) %>% + left_join(tw_info, by=c("grid_id"="id")) +``` + +```{r} +Grid.TW %>% + ggplot(aes(fill = tsai_ratio)) + geom_sf(lwd = 0.1, color="black") + + scale_fill_continuous(high="#2EFF71", low="blue") + + theme_void() + +``` diff --git a/V24_PointMap.Rmd b/V24_PointMap.Rmd new file mode 100644 index 0000000..fb82747 --- /dev/null +++ b/V24_PointMap.Rmd @@ -0,0 +1,112 @@ +```{r include=FALSE} +knitr::opts_chunk$set(cache = T, warning = F, message = F, + class.output = "output", out.width='100%', + fig.asp = 0.5, fig.align = 'center') +library(tidyverse) +``` + +## Mapping Youbike Location + +這個練習的目標是讀取台北市Youbike2.0的站台資料，並繪製點位圖來標記站台滿車或缺車的情形。 + +以下這段程式碼是用於從台北市政府提供的YouBike實時數據API取得資料，並將其轉換為R語言中的資料框架(data frame)格式，最後選擇資料框架的前六個變數來展示。步驟如下： + +1. 載入`httr`和`jsonlite`套件：這兩個套件在R中常用於處理HTTP請求和JSON資料。`httr`套件用於發送網路請求，而`jsonlite`套件則用於解析JSON格式的資料。 + +2. 設定URL：將`url`變數設為指向台北市政府提供的YouBike實時數據API的網址。 + +3. 使用`GET`函數發送請求：透過`httr`套件的`GET`函數向設定的URL發送HTTP GET請求，以取得YouBike的實時數據。 + +4. 解析JSON資料：使用`jsonlite`套件的`fromJSON`函數，將從API獲取到的JSON格式資料解析成R的資料框架。`content`函數用於獲取HTTP回傳的內容，並指定內容格式為"text"，編碼方式為"utf-8"，以確保中文等非ASCII字符能正確顯示。 + +5. 要注意的欄位名稱包含sna（站台名稱）、tot（總車格數）、sbi（現有車數）與lat-lng的經緯度資料。 + +```{r} +library(httr) +library(jsonlite) +url <- "https://tcgbusfs.blob.core.windows.net/dotapp/youbike/v2/youbike_immediate.json" +ubike.df <- fromJSON(content(GET(url),"text", encoding = "utf-8")) +head(ubike.df) %>% select(1:6) + +``` + +### Creating a new variable + +- 為每個站產生一個`sbi`除以`tot`的新變數，來作為ubike站滿不滿的計量單位。 +- 這段程式碼是在前一步的基礎上進一步創建新的 `Fullness` 欄位。該欄位的值根據 `ratio` 欄位的數值判斷，若 `ratio` 大於等於 0.9，則將 `Fullness` 設置為 "Full"；若 `ratio` 小於 0.3，則設置為 "Empty"；否則設置為 "Available"。 + +`if_else` 函數的用法如下： + +``` r +if_else(condition, true_value, false_value) +``` + +- `condition`：一個邏輯向量，用於指定條件判斷。 +- `true_value`：當 `condition` 為真時，返回的值。 +- `false_value`：當 `condition` 為假時，返回的值。 + +`if_else` 函數與 R 原生的 `ifelse` 函數相似，但有一些重要區別。`if_else` 函數在處理缺失值時更加嚴格，要求 `true_value` 和 `false_value` 具有相同的長度和類型，這有助於避免一些潛在的錯誤。 + +```{r} +# ratio <- sbi/tot +ubike.df <- ubike.df %>% + mutate(ratio = sbi/tot) %>% + mutate(Fullness = if_else(ratio >= 0.9, "Full", if_else(ratio < 0.3, "Empty", "Available"))) +``` + +### Mapping with sf + +```{r} +library(ggplot2) +library(sf) + +sf_tpe <- + st_read(dsn = "data/111年9月行政區人口統計_鄉鎮市區_臺北市_SHP/", + layer = "111年9月行政區人口統計_鄉鎮市區", quiet = T) %>% + mutate(across(where(is.character), ~iconv(., from = "BIG5", to = "UTF8"))) %>% + mutate(across(where(is.double), ~if_else(is.na(.),as.double(0),.))) %>% + st_set_crs(3826) %>% st_transform(4326) %>% + filter(COUNTY == "臺北市") + +sf_tpe %>% head() + + +``` + +```{r} +library(sf) +library(ggplot2) +library(dplyr) + +# 你已經有ubike.sf，這是一個sf物件 +ubike.sf <- st_as_sf(ubike.df, coords = c("lng", "lat"), crs = 4326) + +# 使用ggplot繪製地圖和位置點，並根據fullness的值上色 +ggplot() + + geom_sf(data = sf_tpe) + # 繪製sf_tpe區域 + geom_sf(data = ubike.sf, aes(color = Fullness, size=tot), alpha=0.3) + + scale_color_manual(values = c("Full" = "black", "Empty" = "red", "Available" = "blue")) + + theme_void() + +``` + +### Using ggmap (Deprecated) + +- 首先要設定google map的參數，先打開[google map](https://www.google.com.tw/maps)，縮放到你等一下希望看到的地圖底圖範圍後，複製新的網址列如[`https://www.google.com.tw/maps/@25.0353215,121.4916909,12.62z`](https://www.google.com.tw/maps/@25.0353215,121.4916909,12.62z)。這個網址後面的參數包含著經緯度（25.0353215,121.4916909）和地圖的縮放程度（12.62z）。因此你可以把他貼到`get_googlemap()`的函式中。除此之外，尚可以指定地圖種類，自行查詢help（maptype can be terrain, roadmap, satellite, hybrid, or toner-lite）。 + +- 接下來要用`geom_point()`這個`ggplot2`的函式繪製點圖，一共包含四個參數，用`data`指定資料集、用`aes()`指定x、y軸、用`colour`指定顏色、用`size`指定圓圈的大小。這邊我用圓圈的大小來表示總車格數，所以是`ubike.df$tot`，但如果是原始的`tot`數值可能畫起來會太大，這時後就必須要**自行反覆嘗試**，看看要取`log()`或者取`sqrt()`或者除以某個數值，好讓畫出來的點為適合觀看的大小。最後我多加了一個參數`alpha`讓圈圈半透明。 + +```{r eval=FALSE, include=FALSE} +library(ggmap) +library(ggplot2) + +# 設定API金鑰（如果你使用的是Google或Stadia地圖源，需要此步驟） +# register_google(key = "your_google_api_key_here") + +# 獲取地圖數據，這裡以台北市為例，並指定使用OSM作為地圖源 +map_data <- get_map(location = "Taipei", zoom = 12, source = "osm") + +# 使用ggmap來繪製地圖 +ggmap(map_data) + +```