使用 R 語言進行地理計算:入門

2021-02-21 吹鼓手
theme: hpstr highlight: github

knitr::opts_chunk$set(echo = TRUE,
dpi = 400,
warning = FALSE,
message = FALSE,
fig.align = 'center',
comment = "#>")

這是《R 語言地理數據科學》系列課程的第一講，通常這種課程是面向 R 語言高級用戶的，但是考慮到大家的 R 語言可能還沒入門，所以我打算儘可能細緻的講解，讓這門課即使是剛入門的 R 用戶也能輕鬆理解，不過儘管如此，我依舊建議你結合著我的另一個課程：《R 數據科學》一起學習。

本系列的課程的參考資料較多，我還會結合大量實際案例，因此此處我暫時不一一列舉了。我會在每次課程的最後添加參考資料的連結。

為什麼要使用 R 語言進行地理分析？

實際上我並沒有學習過任何桌面 GIS 軟體，例如 QGIS、ArcMap、GRASS、SAGA等，但是程式語言相對桌面軟體的最大優勢就是可重複性了。

所以我也只能從 R 的角度去看這個問題，R 的易用性、跨平臺性以及活躍的社區都是這個問題的答案。基於我自己的使用經驗，使用 R 進行地理計算可以滿足我們 99% 的工作和學習需要。另外結合 R 語言強大的繪圖性能，地理數據可視化也就不是問題了。

例如下圖展示了夜間燈光數據和中國三個城市的位置：

library(leaflet)
popup = c("北京市", "拉薩市", "香港特別行政區")
leaflet(width = "100%", height = "400px") %>%
addProviderTiles("NASAGIBS.ViirsEarthAtNight2012") %>%
addMarkers(lng = c(116.41228, 91.089778, 114.13432),
lat = c(40.185543, 30.036993, 22.381274),
popup = popup)

可重複性：別人可以使用你的代碼和數據得到相同的結果。

R 語言中有諸多用於地理計算和可視化的包，在這裡我主要使用 sf、raster 和 ggplot2 包，另外我會儘可能介紹更多好用的包。我就不再贅述使用 R 語言進行地理分析的歷史了，反正用 sf 包就對了。在數據量不大的時候這個包的性能表現還不錯。

使用 R 語言進行地理分析能做什麼？

這裡展示下我之前的教程中涉及到的一些。

R 中的地理數據準備工作

本節課主要使用下面兩個 R 包：

install.packages('sf')
install.packages('raster')

地理矢量數據

地理矢量數據模型的基礎是位於地理參考系統（crs）中的點。點可以表示為獨立的要素（例如某個公交站臺的位置），也可以組合成複雜的幾何形狀，例如直線和多邊形。

例如 crs = 4326：

library(sf)
library(ggplot2)
library(tidyverse)

worldmp <- read_sf("world.geo.json") %>%
st_transform(4326)

tibble(lon = c(0, 116),
lat = c(0, 40),
name = c("crs：4326 原點", "北京")) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) -> da4326

ggplot(da4326) +
geom_sf_label(aes(label = name), family = cnfont,
nudge_x = -12, nudge_y = -12,
color = "red") +
geom_sf(color = "red", size = 3) +
geom_sf(data = worldmp, fill = NA) +
coord_sf(crs = 4326) +
labs(title = "CRS: 4326") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

我覺得費力的去解釋 4326 意義不大，因為它表示 WGS 84 坐標系，顯然又變成另一個費解的詞語，因為後面我都直接說 crs 的數值了。我國官方要求使用的坐標系是 CCGS 2000，對應的 crs 是：4490，這一坐標系的坐標原點和北京的位置：

worldmp <- read_sf("world.geo.json") %>%
st_transform(4490)

tibble(lon = c(0, 116),
lat = c(0, 40),
name = c("crs：4490 原點", "北京")) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(4490) -> da4490

da4490
da4326
ggplot(da4490) +
geom_sf_label(aes(label = name), family = cnfont,
nudge_x = -12, nudge_y = -12,
color = "red") +
geom_sf(color = "red", size = 3) +
geom_sf(data = worldmp, fill = NA) +
coord_sf(crs = 4490) +
labs(title = "CRS: 4490") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

是不是發現兩個坐標係為什麼差異，所以通常我們不需要特意將 4326 下的坐標轉到 4490 下。兩個坐標系下的坐標僅有細微的差異。

有很多 crs 可供選擇：

# install.packages('rgdal')
rgdal::make_EPSG() %>%
DT::datatable(width = "100%")

sf 包也提供了類似的函數：

sf_proj_info(type = "proj") %>%
DT::datatable(width = "100%")

我們再以 crs = 54003 為例演示下：

worldmp <- read_sf("world_high_resolution_mill.geo.json")
worldmp
# 它的 crs 是 EPSG:54003

rbind(
tibble(lon = 116, lat = 40, name = "北京") %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform(54003),
tibble(X = 0, Y = 0, name = "crs：54003 原點") %>%
st_as_sf(coords = c("X", "Y"), crs = 54003)
) -> da54003

ggplot(da54003) +
geom_sf_text(aes(label = paste0("\n", name)),
family = cnfont,
color = "red") +
geom_sf(color = "red", size = 3) +
geom_sf(data = worldmp, fill = NA) +
coord_sf(crs = 54003) +
labs(title = "CRS: 54003") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

這裡就不再羅列更多的 crs 了，這不是本節課的重點。

simple features

sf 是 simple features 的意思，是代表各種幾何類型的分層數據模型，sf 包支持所有的其中常用幾何類型：

在使用 sf 包進行地理分析的時候，我們通常會將數據讀取為 sf 對象。sf 對象存儲在數據框中，其中地理數據佔據一個特殊列，通常是 geometry，例如我們剛剛讀取的 worldmp:

worldmp
class(worldmp)
plot(worldmp)

通常我使用 ggplot2 繪圖：

ggplot(worldmp) +
geom_sf(aes(fill = name), color = "black",
size = 0.1) +
scale_fill_viridis_d(option = "D") +
theme_ipsum(base_family = cnfont) +
theme(legend.position = "none") +
labs(title = "世界地圖") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

sf 對象的最大優勢就在於，你可以使用幾乎所有處理數據框的手段處理它。例如：

glimpse(worldmp)

為什麼要使用 sfsf 包的函數可以使用 %>% 運算符進行組合，並且可以很好地與 tidyverse 系列的包一起使用。

現在越多越多的 R 包增加了對 sf 對象的支持，例如：tmap 和 mapview。而對於那些還不支持 sf 對象的包，也完全不用擔心，sf 類可以很容易轉換成 sp 類進行使用，例如：

library(sp)
world_sp = as(worldmp, Class = "Spatial")

而 sp 對象也可以很容易的被轉換為 sf 對象：

world_sf = st_as_sf(world_sp)

ggplot2 + sf 的一些基本用法

上面的例子中我們也演示了一些，下面我再通過一個例子演示下：

# 生成一列隨機數
worldmp$pop = runif(215, 1, 100)

# 繪圖
ggplot(worldmp) +
geom_sf(aes(fill = name), color = "black",
size = 0.1) +
stat_sf_coordinates(geom = "point", color = "#E31A1C",
aes(size = pop), shape = 1) +
scale_fill_viridis_d(option = "D") +
theme_ipsum(base_family = cnfont) +
theme(legend.position = "none") +
labs(title = "世界地圖") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

幾何類型：簡單要素幾何（sfg）

geometry 是 sf 的基礎，本講中我們將一一介紹常見的七種類：POINT，LINESTRING，POLYGON，MULTIPOINT，MULTILINESTRING，MULTIPOLYGON 和 GEOMETRYCOLLECTION。

POINT、LINESTRING 和 POLYGON

# POINT：點
st_point(c(5, 2))

# LINESTRING：線
st_linestring(matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T))

# POLYGON：多邊形
st_polygon(list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
ncol = 2, byrow = T)))

MULTI～

# MULTIPOINT
st_multipoint(matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T))

# MULTILINESTRING
st_multilinestring(list(
matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T),
matrix(c(1, 2, 2, 4),
ncol = 2, byrow = T)
))

# MULTIPOLYGON
st_multipolygon(list(
list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
ncol = 2, byrow = T)),
list(matrix(c(0, 2, 1, 2, 1, 3, 0, 3, 0, 2),
ncol = 2, byrow = T))
))

GEOMETRYCOLLECTION: 大雜燴

st_geometrycollection(list(
st_point(c(5, 2)),
st_linestring(matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T)),
st_polygon(list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
ncol = 2, byrow = T))),
st_multipoint(matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T)),
st_multilinestring(list(
matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
ncol = 2, byrow = T),
matrix(c(1, 2, 2, 4),
ncol = 2, byrow = T)
)),
st_multipolygon(list(
list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
ncol = 2, byrow = T)),
list(matrix(c(0, 2, 1, 2, 1, 3, 0, 3, 0, 2),
ncol = 2, byrow = T))
))
)) %>% plot()

幾何類型：簡單要素列表（sfc）

point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc
# 還可以設定 crs
st_sfc(point1, point2, crs = 4326)

sf

我們可以按照 sfg -> sfc -> sf 的順序構造 sf 對象，這樣的操作在實際中很少使用，僅僅用於演示 sf 對象是如何構建的：

# sfg
bj_point <- st_point(c(116, 40))
# sfc
bj_geom <- st_sfc(bj_point, crs = 4326)
# tibble
bj <- tibble(name = "北京")
bj_sf = st_sf(bj, geometry = bj_geom)
bj_sf

# 繪圖展示
ggplot(worldmp) +
geom_sf(aes(fill = name), color = "black",
size = 0.1) +
geom_sf(data = bj_sf, color = "red") +
scale_fill_viridis_d(option = "D") +
theme_ipsum(base_family = cnfont) +
theme(legend.position = "none") +
labs(title = "世界地圖") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())

raster：柵格數據

地理柵格數據通常由定義坐標參考系統，範圍和原點的 header 和代表像素的矩陣組成。我平時遇到 raster 數據都是簡單處理之後轉為 sf 對象再進行處理，例如夜間燈光數據（2013年全球的），下面的示例很費時間，建議在夜間運行：

F182013.v4c_web.avg_vis.tif 數據下載自：https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html ，為了方便大家下載，我提供一個百度網盤連結：
連結:https://pan.baidu.com/s/16RmU9FHtFbO5i4aWb_vujg 密碼:b2cv

library(raster)
raster('F182013.v4c_web.avg_vis.tif') -> light

# 從其中提取中國的
worldmp %>% 
  dplyr::filter(name == "中國") %>% 
  st_transform(projection(light)) -> cn
mask(light, cn) -> cnlight

# 將 raster 數據轉為 sf 
cnlight %>% 
  rasterToPoints(spatial = TRUE) %>% 
  st_as_sf() -> cnlight_sf

# 繪圖展示
ggplot(cnlight_sf) + 
  geom_sf(aes(color = F182013.v4c_web.avg_vis), size = 0.01) + 
  geom_sf(data = cn, size = 0.3, color = "white", fill = NA) + 
  scale_color_gradientn(
    colours = c("#1E1E1E", "#1E1E1E", "#fcfff7", "#f1ff8f", "#f1e93f", "#fabd08")) + 
  theme_modern_rc(base_family = cnfont,
                  plot_title_family = cnfont,
                  subtitle_family = cnfont,
                  caption_family = cnfont) + 
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) + 
  labs(title = "中國夜間燈光數據",
       subtitle = "RStata Project",
       caption = "數據來源：Earth Observation Group | ngdc.noaa.gov\n<https://www.ngdc.noaa.gov/eog/dmsp/downloadV4composites.html>")

坐標參考系

地理坐標系使用兩個值（經度和緯度）識別地球表面上的任何位置。投影坐標參考系統則是將地球表面展開得到的笛卡爾坐標系，不同的展開方式就得到了不同投影下的坐標系。投影的缺點就是必然產生失真（面積、距離、形狀和方向）。

R 中描述 crs 的方法主要是使用 epsg 代碼和 proj4string。這兩種方法都有優點和缺點。epsg 代碼通常較短，更容易記住，但是該代碼僅引用一種定義明確的坐標參考系；proj4string 則可以更靈活的定義各種參數（例如投影類型，基準面和橢球體）。

我們可以使用 st_crs 查看某個 sf 對象的 crs：

st_crs(worldmp)

# 查看 proj4string
st_crs(worldmp)$proj4string

# 查看度量單位
st_crs(worldmp)$units

# 查看 EPSG 代碼
st_crs(worldmp)$epsg

如果某個對象缺少 crs 或者設置了錯誤的 crs，你可以使用 st_set_crs 設定

st_set_crs(worldmp, 54003)

projection 函數可以提取 raster 對象的 crs：

projection(light)
#> [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

單位

我們先了解一下 R 語言裡面的測度：

library(units)

長度

l <- set_units(1:10, m)
l

我們再把米轉換成釐米：

set_units(l, cm)

還可以進行計算：

set_units(l, cm) + l

面積

例如 400 公頃等於 4 平方千米：

a <- set_units(400, ha)
# 等於 4 平方千米
set_units(a, km2)

矢量

例如速度的定義和轉換：

vel <- set_units(seq(20, 50, 10), km/h)
set_units(vel, m/s)

sf 對象的一個優勢就是對單位的支持，例如我們計算世界各國的面積：

st_area(worldmp)

我們可以把它合併到 worldmp 數據框裡面：

worldmp$area <- set_units(st_area(worldmp), km2)
worldmp

作業：

如果時間允許的話，可以試著從全球的夜間燈光數據裡提取美國的部分並繪圖展示。參考資料Chapter 2 Geographic data in R | Geocomputation with R. https://geocompr.robinlovelace.net/spatial-class.html

相關焦點

  • R 語言入門學習路線與資源匯總
    剛開始學習R語言，可以看一些免費入門視頻，可以在B站上面搜，B站有很多免費的R語言入門視頻，比如《尚學堂尹鴻的R語言速成實戰》[1]，這個視頻是免費的，視頻講義素材來源於《R語言實戰》的第一章到第六章，這幾章可以邊看書邊跟著視頻學習，用來簡單入門基本夠了。
  • R 語言入門精講 Day#1
    這種程式語言被命名為R語言，基於兩個R語言作者的名字的第一個字母（Robert Gentleman和Ross Ihaka），並且部分是貝爾實驗室語言S的名稱。[1]簡單來說，R是一門統計計算語言，是一個用於統計計算和統計製圖的優秀工具，是一套開源的數據分析解決方案。R 語言的特點如前所述，R語言是用於統計分析，圖形表示和報告的程式語言和軟體環境。
  • 想自學R語言,從哪裡入門|R教程
    最令人沮喪的是——許多統計分析無法實現，或滯後，例如結構方程模型、meta分析（元分析）等。這意味著為了實現不同的功能，我要在內存有限的情況下，下載1+N個軟體。可能一不小心，還要耗費更多經濟成本（購買軟體）或時間成本（破解版）。綜合可視化 + 統計需求兩個原因，我終於對R下手了。雖然學習的時間並不長，大概只是脫離了小白的階段，但還是厚著臉皮來做個分享。
  • 地理探測器R語言包使用方法(中文版)
    地理探測器「GD」R包使用說明（中文）Song, Y., Wang, J., Ge, Y.關於「GD」包地理探測器可解決以下問題：從空間異質性角度探索潛在因素或解釋變量。探索地理變量的潛在交互影響。從潛在的解釋變量中識別高風險或低風險區域。
  • 為什麼要使用R語言?我的R語言學習盤點
    R可以像MATLAB一樣解決有關數值計算的問題，而且具有強大的數據處理，繪圖功能。R擁有大量的統計分析工具包，我的感覺是——只有我們沒聽說過的工具，絕對沒有R沒有的工具包。配合著各種各樣的工具包，你可以毀滅任何關於數據和統計的問題。因為數據包的數量龐大，所以查找自己需要的數據包，可能很煩惱。
  • R語言入門
    後臺有一些朋友希望我能多介紹一下R語言在金融中的應用，這篇我就基於自身經驗簡單談談如何入門R語言吧。
  • R語言之控制流的使用方法
    關於這兩種分支我們看下圖：對於二分支，我個人更喜歡使用ifelse()函數,該函數是R語言內置的函數，運行速度是非常快的。雖然使用within()函數實現重編碼有一點點複雜，但計算效率還是非常高的！在知道循環多少次的前提下可使用這樣的for語句。上面這個例子不是一對一的重編碼問題，建議使用within()內置函數解決編碼，可以大大提高計算速度。
  • R語言入門課程
    針對R語言小白的入門課程，已經進行了半個多月：R語言小白的免費入門課程。現在對講述的內容簡單做個梳理。
  • R語言學習入門導航-特別版
    r-bloggers: How to learn Rr-bloggers: Tutorials for learning RR Learning Path: From beginner to expert in R in 7 stepsRstudio online-learningBookdown Homew3cschool
  • 【R語言】入門知識--常用操作和例子
    2 簡單的數學計算、數學函數以及如何編輯R程序 1，用R進行簡單的計算： 我們可以用R進行以下各類運算。例如： 輸入： sqrt(2) 結果：[1] 1.414213 你還可以使用多個函數進行計算。例如： 輸入： sqrt(100) + round(100) / log10(100) 結果：[1] 60 下面介紹的是R中可以使用的數學函數。
  • Tidy時代R語言學習的一些ABC
    作為一種統計程式語言，R最初的用戶主要集中在學術領域，開源的優勢和快速、便利的統計建模，可以讓研究人員擺脫昂貴又不夠靈活的各種商業軟體。大數據革命的出現，將R帶到了數據科學舞臺的中央，R語言從學術界破圈而出，獲得了迅猛的發展。然而，R本身作為統計語言的一些固有缺點在越來越複雜的數據問題面前也逐漸暴露出來。比如，R的語法不夠規範，計算效率慢、工業化部署能力差等等。
  • R語言安裝及入門
    當然，若是一款軟體既能統計，又能夠繪圖，那肯定是再好不過的啦~簡單來說，R就是這樣一門擅長於統計，數據可視化尤其強大的語言。它既可以統計，也可以進行繪圖，最重要的是它是完全開源和免費的！那今天就請Ten years old詳細講解一下，如何進行R的安裝及入門？
  • R語言 | Tidyverse包入門介紹
    本期將介紹以下函數運算：%>%, tibble, as_tibble, read_csv, read_delim, filter, select, mutate, group_by, summarise, gather, spread數據處理流程 (R for Data Science, pp. ix)迄今為止，我們的R語言入門級講解，都會對已有的數據進行預處理
  • 地理數據科學培訓班第一課之初識R語言
    歡迎收看明晚九點的直播講解：R 和 RStudio 的安裝及 R Profile 的配置 & 初識 R 語言數據爬取本次直播講解將使用 Windows 系統的電腦演示。，下載這個安裝：https://cran.r-project.org/bin/windows/Rtools/rtools40-i686.exe安裝完成之後還需要進行環境變量的配置，稍後我再介紹。
  • 如何在 Pycharm 中高效使用 R 語言 (圖文詳解)
    相信大家學習生信的時候，都會或多或少使用
  • C語言小白入門第一課
    C 語言是一種通用的、面向過程式的電腦程式設計語言。1972 年，為了移植與開發 UNIX 作業系統，丹尼斯·裡奇在貝爾電話實驗室設計開發了 C 語言。C 語言是一種廣泛使用的計算機語言，它與 Java 程式語言一樣普及，二者在現代軟體程式設計師之間都得到廣泛使用。
  • R語言的三種聚類方法
    一、層次聚類1)距離和相似係數r語言中使用dist(x, method =
  • R語言包,掌握這10個就夠了!(含資源下載)
    ggplot2做出來的圖更漂亮，但它入門難啊~lattice則適合入門選手，作圖速度較快，還能進行三維繪圖，這是ggplot2不具備的。zooforecast裡有很多好用的內置函數，但有時候還需要更簡單易用的計算移動平均、移動標準差的函數，這時zoo會是個很好的補充。spatstat空間分析神器，比如某種疾病發病率在地理空間上的分布特徵等等。
  • R語言中使用subset函數對數據進行分類管理操作
    我們在SCI論文中常常可以見到這樣的表格，是根據分類來做出統計結果的，如下圖，是根據患者是否存活把患者分成了兩類倖存的和死亡的做分別統計，然後得出各類統計結果那麼，R語言是怎麼做出這樣的表格呢首先我們要把數據進行分割，得到一個倖存的數據表和一個死亡的數據表，然後再分別統計，我們今天利用R語言自帶的subset函數來演示這一功能，這是一個非常重要的功能，為今後我們對數據進一步分析做準備。
  • R語言向量化運算:apply函數族用法心得
    R語言和Python的忠實擁躉，為成為一名未來的數據科學家而奮鬥終生。個人公眾號：數據科學家養成記 （微信ID：louwill12）當初入坑R語言的時候，就在各種場合看到老司機的忠告，「儘量避免使用循環！」一開始並不明白這其中的奧義，直到後來對R語言有深入接觸後，才領會R語言在向量化運算方面的強大功能。本篇內容就總結小編在使用R語言向量化運算apply函數族的一些心得體會。