作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具(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)

参考资料：