作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具（R），致力于利用数据科学迅速积累行业经验优势和科学知识发现，涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等，著有《R语言高效数据处理指南》（《R语言数据高效处理指南》(黄天元)
本帖子会简单介绍如何读入并处理栅格数据。首先，我们会用到一个矢量数据，数据来自：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/"
)
# 读取特定图层
Ausoutline <- st_read("data/",
layer='gadm36_AUS_0')
可以对这个数据集进行观察：
# 观察
print(Ausoutline)
# 查看proj4坐标系
st_crs(Ausoutline)$proj4string
然后，让我们读入栅格数据，这是一个全球平均气温数据，来源于：Historical climate data（tavg 5m）。
# 读入栅格数据
jan<-raster( "data/")
# have a look at the raster layer jan
jan
plot(jan)
我们可以改变它的坐标系，然后再次绘图：
# set the proj 4 to a new object
newproj<-"+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# get the jan raster and give it the new proj4
pr1 <- jan %>%
projectRaster(., crs=newproj)
plot(pr1)
下面，让我们一次性读入所有栅格数据：
# 一次性读入所有栅格数据
dir("data/",full.names = T) %>%
stack() -> worldclimtemp
worldclimtemp
让我们为每一个图层命名，它其实是每个月的平均温度，因此可以用月份命名：
month <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
names(worldclimtemp) <- month
如果想取出一月份的数据，有两种方法：
worldclimtemp[[1]]
worldclimtemp$Jan
下面，让我们提取特定城市的数据，这是通过经纬度来提取的：
## 提取特定城市的数据
site <- c("Brisbane", "Melbourne", "Perth", "Sydney", "Broome", "Darwin", "Orange",
"Bunbury", "Cairns", "Adelaide", "Gold Coast", "Canberra", "Newcastle",
"Wollongong", "Logan City" )
lon <- c(153.03, 144.96, 115.86, 151.21, 122.23, 130.84, 149.10, 115.64, 145.77,
138.6, 153.43, 149.13, 151.78, 150.89, 153.12)
lat <- c(-27.47, -37.91, -31.95, -33.87, 17.96, -12.46, -33.28, -33.33, -16.92,
-34.93, -28, -35.28, -32.93, -34.42, -27.64)
#Put all of this inforamtion into one list
samples <- data.frame(site, lon, lat, "site")
samples
# Extract the data from the Rasterstack for all points
AUcitytemp<- raster::extract(worldclimtemp, samples)
提取的变量是一个矩阵，可以稍微进行转化，把城市名称也加上。
AUcitytemp
Aucitytemp2 <- AUcitytemp %>%
as_tibble()%>%
add_column(Site = site, .before = "Jan")
Aucitytemp2
现在，我们要从世界气温图中取出澳洲的部分，实现如下（利用crop函数）：
Austemp <- Ausoutline %>%
# now crop our temp data to the extent
crop(worldclimtemp,.)
# plot the output
plot(Austemp)
这个图中还包括了非澳洲的部分，如果只想要澳洲的部分，可以这样操作：
exactAus <- Austemp %>%
mask(Ausoutline, )
plot(exactAus)
可以尝试对三月份的数据做一个温度分布直方图：
exactAusdf <- exactAus %>%
()
# set up the basic histogram
gghist <- ggplot(exactAusdf,
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,
)),
color="blue",
linetype="dashed",
size=1)+
theme(plot.title = element_text(hjust = ))
更多可视化的选择可参考：https://andrewmaclachlan.github.io/CASA0005repo/rasters-descriptive-statistics-and-interpolation.html#histogram-with-ggplot
最后，我们会演示一下空间插值的过程。首先，我们来合并气温的时空分布数据：
samplestemp<-AUcitytemp%>%
cbind(samples)
samplestemp
而后，让我们把它转换为空间信息数据框，这里需要声明经纬所在列，以及坐标系（这里定义为4326，也就是WGS 84）：
samplestemp<-samplestemp%>%
st_as_sf(.,coords = c("lon", "lat"),
crs =4326,
agr = "constant")
samplestemp
这里我们可以对澳大利亚轮廓图和城市分布进行可视化：
plot(Ausoutline$geom)
plot(st_geometry(samplestemp), add=TRUE)
插值之前，让我们先统一坐标系，然后转换为SP对象：
samplestemp <- samplestemp %>%
st_transform(3112)
Ausoutline <- Ausoutline %>%
st_transform(3112)
samplestempSP <- samplestemp %>%
as(., 'Spatial')
AusoutlineSP <- Ausoutline %>%
as(., 'Spatial')
现在意味着，我们要用手头若干个城市的数值，来对整个澳大利亚的气温空间分布进行插值。我们需要创建一个要插值的空间范围：
emptygrd <- (spsample(AusoutlineSP, n=1000, type="regular", cellsize=200000))
names(emptygrd) <- c("X", "Y")
coordinates(emptygrd) <- c("X", "Y")
gridded(emptygrd) <- TRUE # Create SpatialPixel object
fullgrid(emptygrd) <- TRUE # Create SpatialGrid object
# Add the projection to the grid
proj4string(emptygrd) <- proj4string(samplestempSP)
然后进行插值：
p_load(gstat)
# Interpolate the grid cells using a power value of 2
interpolate <- gstat::idw(Jan ~ 1, samplestempSP, newdata=emptygrd, idp=)
这里用的IDW算法来插值，只使用1月份的数据。idp参数越大，则随着距离的增大，数值减少越大。如果idp为0，那么随着距离加大，依然不会有任何数值衰减。关于方法细节，可参考：How inverse distance weighted interpolation works。最后，我们看一下插值之后的结果：
# Convert output to raster object
ras <- raster(interpolate)
# Clip the raster to Australia outline
rasmask <- mask(ras, Ausoutline)
# Plot the raster
plot(rasmask)
