作者：黃天元，復旦大學博士在讀，熱愛數據科學與開源工具(R)，致力于利用數據科學迅速積累行業經驗優勢和科學知識發現，涉獵內容包括但不限于信息計量、機器學習、數據可視化、應用統計建模、知識圖譜等，著有《R語言高效數據處理指南》(《R語言數據高效處理指南》(黃天元)【摘要 書評 試讀】- 京東圖書)。知乎專欄：R語言數據挖掘。郵箱：huang.tian-yuan@qq.com.歡迎合作交流。HopeR：R語言空間數據分析(零)：總目錄?zhuanlan.zhihu.com

本帖子會簡單介紹如何讀入并處理柵格數據。首先，我們會用到一個矢量數據，數據來自：https://gadm.org/download_country_v3.html，用到的是澳洲的地圖。讀取方法如下：

# 獲得數據的方法之一

# wget --no-check-certificate https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_AUS_gpkg.zip

library(pacman)

p_load(sf,raster,tidyverse)

# 查看有哪些圖層

st_layers(

"data/gadm36_AUS.gpkg"

)

# 讀取特定圖層

Ausoutline

layer='gadm36_AUS_0')

可以對這個數據集進行觀察：

# 觀察

print(Ausoutline)

# 查看proj4坐標系

st_crs(Ausoutline)$proj4string

plot(Ausoutline)

然后，讓我們讀入柵格數據，這是一個全球平均氣溫數據，來源于：Historical climate data(tavg 5m)。

# 讀入柵格數據

jan

# have a look at the raster layer jan

jan

plot(jan)

我們可以改變它的坐標系，然后再次繪圖：

# set the proj 4 to a new object

newproj

# get the jan raster and give it the new proj4

pr1 %

projectRaster(., crs=newproj)

plot(pr1)

下面，讓我們一次性讀入所有柵格數據：

# 一次性讀入所有柵格數據

dir("data/wc2.1_5m_tavg",full.names = T) %>%

stack() -> worldclimtemp

worldclimtemp

讓我們為每一個圖層命名，它其實是每個月的平均溫度，因此可以用月份命名：

month

"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

names(worldclimtemp)

如果想取出一月份的數據，有兩種方法：

worldclimtemp[[1]]

worldclimtemp$Jan

下面，讓我們提取特定城市的數據，這是通過經緯度來提取的：

## 提取特定城市的數據

site

"Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle",

"Wollongong", "Logan City" )

lon

138.6, 153.43, 149.13, 151.78, 150.89, 153.12)

lat

-34.93, -28, -35.28, -32.93, -34.42, -27.64)

#Put all of this inforamtion into one list

samples

samples

# Extract the data from the Rasterstack for all points

AUcitytemp" />

提取的變量是一個矩陣，可以稍微進行轉化，把城市名稱也加上。

AUcitytemp

Aucitytemp2 %

as_tibble()%>%

add_column(Site = site, .before = "Jan")

Aucitytemp2

現在，我們要從世界氣溫圖中取出澳洲的部分，實現如下(利用crop函數)：

Austemp %

# now crop our temp data to the extent

crop(worldclimtemp,.)

# plot the output

plot(Austemp)

這個圖中還包括了非澳洲的部分，如果只想要澳洲的部分，可以這樣操作：

exactAus %

mask(Ausoutline, na.rm=TRUE)

plot(exactAus)

可以嘗試對三月份的數據做一個溫度分布直方圖：

exactAusdf %

as.data.frame()

# set up the basic histogram

gghist

aes(x=Mar)) +

geom_histogram(color="black",

fill="white")+

labs(title="Ggplot2 histogram of Australian March temperatures",

x="Temperature",

y="Frequency")

# add a vertical line to the hisogram showing mean tempearture

gghist + geom_vline(aes(xintercept=mean(Mar,

na.rm=TRUE)),

color="blue",

linetype="dashed",

size=1)+

theme(plot.title = element_text(hjust = 0.5))

最后，我們會演示一下空間插值的過程。首先，我們來合并氣溫的時空分布數據：

samplestemp%

cbind(samples)

samplestemp

而后，讓我們把它轉換為空間信息數據框，這里需要聲明經緯所在列，以及坐標系(這里定義為4326，也就是WGS 84)：

samplestemp%

st_as_sf(.,coords = c("lon", "lat"),

crs =4326,

agr = "constant")

samplestemp

這里我們可以對澳大利亞輪廓圖和城市分布進行可視化：

plot(Ausoutline$geom)

plot(st_geometry(samplestemp), add=TRUE)

插值之前，讓我們先統一坐標系，然后轉換為SP對象：

samplestemp %

st_transform(3112)

Ausoutline %

st_transform(3112)

samplestempSP %

as(., 'Spatial')

AusoutlineSP %

as(., 'Spatial')

現在意味著，我們要用手頭若干個城市的數值，來對整個澳大利亞的氣溫空間分布進行插值。我們需要創建一個要插值的空間范圍：

emptygrd

names(emptygrd)

coordinates(emptygrd)

gridded(emptygrd)

fullgrid(emptygrd)

# Add the projection to the grid

proj4string(emptygrd)

然后進行插值：

p_load(gstat)

# Interpolate the grid cells using a power value of 2

interpolate

這里用的IDW算法來插值，只使用1月份的數據。idp參數越大，則隨著距離的增大，數值減少越大。如果idp為0，那么隨著距離加大，依然不會有任何數值衰減。關于方法細節，可參考：How inverse distance weighted interpolation works。最后，我們看一下插值之后的結果：

# Convert output to raster object

ras

# Clip the raster to Australia outline

rasmask

# Plot the raster

plot(rasmask)

參考資料：

總結

以上是生活随笔為你收集整理的r语言worldclim数据_R语言空间数据分析（五）：栅格数据处理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。