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

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

# 獲得數(shù)據(jù)的方法之一

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

可以對(duì)這個(gè)數(shù)據(jù)集進(jìn)行觀察：

# 觀察

print(Ausoutline)

# 查看proj4坐標(biāo)系

st_crs(Ausoutline)$proj4string

plot(Ausoutline)

然后，讓我們讀入柵格數(shù)據(jù)，這是一個(gè)全球平均氣溫?cái)?shù)據(jù)，來源于：Historical climate data(tavg 5m)。

# 讀入柵格數(shù)據(jù)

jan

# have a look at the raster layer jan

jan

plot(jan)

我們可以改變它的坐標(biāo)系，然后再次繪圖：

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

下面，讓我們一次性讀入所有柵格數(shù)據(jù)：

# 一次性讀入所有柵格數(shù)據(jù)

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

stack() -> worldclimtemp

worldclimtemp

讓我們?yōu)槊恳粋€(gè)圖層命名，它其實(shí)是每個(gè)月的平均溫度，因此可以用月份命名：

month

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

names(worldclimtemp)

如果想取出一月份的數(shù)據(jù)，有兩種方法：

worldclimtemp[[1]]

worldclimtemp$Jan

下面，讓我們提取特定城市的數(shù)據(jù)，這是通過經(jīng)緯度來提取的：

## 提取特定城市的數(shù)據(jù)

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

提取的變量是一個(gè)矩陣，可以稍微進(jìn)行轉(zhuǎn)化，把城市名稱也加上。

AUcitytemp

Aucitytemp2 %

as_tibble()%>%

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

Aucitytemp2

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

Austemp %

# now crop our temp data to the extent

crop(worldclimtemp,.)

# plot the output

plot(Austemp)

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

exactAus %

mask(Ausoutline, na.rm=TRUE)

plot(exactAus)

可以嘗試對(duì)三月份的數(shù)據(jù)做一個(gè)溫度分布直方圖：

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

最后，我們會(huì)演示一下空間插值的過程。首先，我們來合并氣溫的時(shí)空分布數(shù)據(jù)：

samplestemp%

cbind(samples)

samplestemp

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

samplestemp%

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

crs =4326,

agr = "constant")

samplestemp

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

plot(Ausoutline$geom)

plot(st_geometry(samplestemp), add=TRUE)

插值之前，讓我們先統(tǒng)一坐標(biāo)系，然后轉(zhuǎn)換為SP對(duì)象：

samplestemp %

st_transform(3112)

Ausoutline %

st_transform(3112)

samplestempSP %

as(., 'Spatial')

AusoutlineSP %

as(., 'Spatial')

現(xiàn)在意味著，我們要用手頭若干個(gè)城市的數(shù)值，來對(duì)整個(gè)澳大利亞的氣溫空間分布進(jìn)行插值。我們需要?jiǎng)?chuàng)建一個(gè)要插值的空間范圍：

emptygrd

names(emptygrd)

coordinates(emptygrd)

gridded(emptygrd)

fullgrid(emptygrd)

# Add the projection to the grid

proj4string(emptygrd)

然后進(jìn)行插值：

p_load(gstat)

# Interpolate the grid cells using a power value of 2

interpolate

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

# Convert output to raster object

ras

# Clip the raster to Australia outline

rasmask

# Plot the raster

plot(rasmask)

參考資料：

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。