前言

R语言作为一门开源的统计分析及可视化语言，目前已经拥有许多涵盖众多领域的包，其中就包括在空间数据领域里很出门的sf包和Raster包，它们能够很好地处理矢量和栅格两种数据，那么这次我们就一起来学习一下运用这些包来处理空间数据。
当然我也是站在前人的肩膀上学习的，因此这次的资料很多地方也是借鉴过来的~

知识目录

• 栅格数据读取
• 栅格数据可视化
• 栅格数据的处理
  • 计算平均值、最大值、最小值、总和
  • 栅格重分类
  • 栅格裁剪
  • 栅格掩膜
• 矢量数据读取
• 矢量数据可视化
• 地图制图
• 空间数据的保存

栅格数据读取

数据介绍
本次实验的数据是来自中科院地理所的逐年降水量和逐年平均气温数据，单位分别是0.1毫米和0.1摄氏度。
这里我选取了连续的三个年份的TIF格式数据来实验，方法是不变的。


加载相应的库

library(raster) #处理栅格数据
library(rasterVis) # 栅格数据可视化制图
library(ggspatial) # ggplot2的空间数据拓展包，用于可视化
library(ggplot2)   # 最经典的R语言绘图包之一
library(sf) #读取矢量数据常用的包

使用raster()函数读取tiff文件
这里我的数据存放在了data的文件夹里，因此Raster_Path[2]的意思是该路径下的第2个数据，即1990年的气温数据

#设置Tif数据存放的路径
Raster_Path <- dir("D:/New Desktop/Temp/data",full.names = TRUE)
#使用raster()函数读取路径下指定序号的栅格数据
Temp1990 <- raster(Raster_Path[2])
#使用plot()绘图函数进行可视化
plot(Temp1990)
#使用levelplot()函数进行可视化，levelplot()函数来自于rasterVis包
levelplot(Temp1990, #绘图对象
          margin = list(FUN = 'median'), #设置边缘分布
          par.settings=RdBuTheme) #颜色主题设置

多个栅格数据可视化
使用stack()函数可以将路径下指定序号的文件合成一个栅格数据集

#读取多个tif
Temp_all <- stack(Raster_Path[c(1,2,3)])
plot(Temp_all) # 使用基础的plot函数
levelplot(Temp_all)  # 使用levelplot函数

栅格数据处理

1、栅格数据的四则运算

这里分为两种情况，一种情况是单份栅格数据的四则运算，第二种情况是多份栅格数据的四则运算，后者往往是我们处理长时序数据需要用到的操作
单份栅格数据运算——calc()函数
由于本文中的实验气温数据单位是0.1摄氏度，所以我这里就将气温数据除以10，修正放大后的摄氏度。

calc(x,fun)函数
x是需要进行运算的栅格数据
fun是自定义的函数

#将第一个年份的栅格值缩小10倍
Temp1990_correct <- calc(Temp1990,function(x) x/10)
plot(Temp1990_correct)

多份栅格数据运算

一、使用overlay()函数
overlay(x, y, …, fun, filename = “”, recycle = TRUE, forcefun = FALSE)
x、y、…：参与叠加的栅格对象；
fun：进行属性数据运算的函数。
二、简单粗暴的直接运算
直接把栅格数据作为操作对象进行加减乘除

Temp1990 <- raster(Raster_Path[1])
Temp1991 <- raster(Raster_Path[2])
Temp1992<- raster(Raster_Path[3])
#使用overlay()函数
Temp_sub <- overlay(Temp1990,Temp1991, 
                    fun = function(x, y) (x+y)/2) #求两个年份的气温平均值
#直接进行运算
result <- (Temp1990+Temp1991)/3
plot(Temp_sub)
plot(result)

2、栅格数据的均值、最值、求和计算

• mean() 平均值
• max() 最大值
• min() 最小值
• sum() 求和
• cellStats(x,stat=‘mean’) 栅格统计
  x表示待计算的栅格；
  stat表示需要统计的内容，包括sum, mean, min, max, sd（标准差）, skew（偏度） 和 rms（均方根）；

#使用mean()函数求取平均值，也可以换成max、sum、sd等等
Temp_ave <- mean(Temp_all,na.rm=TRUE)
plot(Temp_ave)

3、栅格重分类

栅格的重分类需要用到一份重分类矩阵以及reclassfy()函数,重分类矩阵需要定义栅格数值的区间和重分类后的编号，类似于“start”，“end”，“ID”这样，具体看以下例子

#第一列表示起始数值，第二列表示终止数值，第三列的1、2、3、4是重分类的编号
# 建立重分类矩阵
reclassmatrix <- matrix(c(-159,-150,1,
                          -150,-100,2,
                          -100,-50,3,
                          -50,0,4,
                          0,50,5,
                          50,100,6,
                          100,150,7,
                          150,200,8,
                          200,250,9,
                          250,270,10
                          ),10,3,byrow = TRUE)
Temp_reclass <- reclassify(Temp1990,rcl = reclassmatrix)
plot(Temp_reclass)

好吧，我承认这个重分类的效果不是很好，估计是我重分类的矩阵设置得不是很好，但方法应该没有问题。

矢量数据的读取

这里为什么突然讲到了矢量数据的读取呢？主要是因为后面的裁剪和掩膜需要用到矢量数据，所以这里简单介绍一下读取和可视化矢量数据的方法。
读取矢量数据主要用到sf包里的read_sf()函数，用法如下所示：

# 读取裁剪范围的矢量数据
area <- read_sf("D:/New Desktop/GISdata/粤港澳大湾区.shp")
plot(area)

那看完矢量数据的可视化之后，有小伙伴就会问怎么它把每个字段都绘制了一遍，难道不能只绘制出它的形状吗？
这个问题其实跟sf包存储空间数据的格式有很大关系，因为它是将空间数据的几何信息和属性信息分开存储的，几何信息存储在geometry里，属性信息则是以一列列的字段形式存储，所以如果我们只想浏览矢量数据的形状，那么可以用代码plot(data$geometry)来可视化它的几何信息。

#浏览矢量数据的几何信息
plot(area$geometry)

5、栅格裁剪(Crop)

栅格裁剪需要用一个矢量数据或者栅格栅格对待操作的数据进行裁剪。
用到的主要是crop()函数

切片函数crop()使用矢量数据的矩形范围对栅格数据进行裁剪。函数语法结构如下：
crop(x, y, filename = “”, snap = ‘near’, datatype = NULL, …)

• x：待裁剪的栅格对象；
• y：裁剪范围或包含裁剪范围的空间数据对象（矢量、栅格数据均可）；
• filename：裁剪后的数据储存地址，可选；
• snap：数据对齐方式，可选项有near、in、out。

# 读取裁剪范围的矢量数据
area <- read_sf("D:/New Desktop/GISdata/粤港澳大湾区.shp")
plot(area)
#查看矢量要素的投影信息
crs(area)
#查看栅格数据的投影信息
crs(Temp1990)
#投影信息转换
area_wgs <- st_transform(area,4326)
#查看转换后的投影信息
crs(area_wgs)
clip <- crop(Temp1990,area_wgs,snap = "near")
plot(clip)

6、栅格掩膜(Mask)

掩膜函数mask() 保持矢量数据实际范围内的属性数据不变，而将范围外的属性数据都转为其他值（默认转换为空值NA）。函数语法结构如下：

mask(x, mask, filename = “”, inverse = FALSE, maskvalue = NA, updatevalue = NA,updateNA = FALSE, …)

• x：栅格对象；
• mask：作为掩膜的空间数据对象（矢量、栅格数据均可）；
• inverse：是否执行反向掩膜操作；
• maskvalue：掩膜内的栅格单元的属性值若为该值则将其转换为updatevalue参数对应的数值；
• updatevalue：掩膜外的所有栅格单元的属性值（空值除外）会被转换为该值，默认为空值NA；
• updateNA：若为TRUE，则掩膜外的空值也会被转换为updatevalue参数对应的数值。

mask <- mask(Temp1990,area_wgs,inverse = FALSE)
plot(mask,main = "Mask Data")#顺手添加个标题

裁剪与掩膜需要注意的地方

1. 切片会改变原有栅格对象的范围，但其得到的是一个矩形范围
2. 掩膜不会改变原有栅格对象的范围，仅是将矢量范围外的数据变为空值，因此右图看起来大湾区的范围很小
3. 对切片后的栅格对象再进行掩膜操作，也可以得到预期裁剪效果

下面是对裁剪完之后的栅格数据进行掩膜提取，得到当前比例尺下所需的区域

mask_area <- mask(clip,area_wgs)
plot(mask_area)

关于数据投影信息的问题

在上面的栅格数据操作过程中，需要注意的地方是进行操作的数据需要保持范围一致，分辨率一致等，因此我们时常需要检查加载进来的数据的投影信息，必要时需要进行坐标系转换
投影转换的方法分为矢量数据和栅格数据两种:
对于矢量数据：

data_proj <- st_transform(data,crs(data2))    #将data转换为data2的投影，并保存为data_proj
data_wgs <- st_taranform(data,4326)  #将数据转换为WGS84地理坐标系

对于栅格数据：

RasterData <- projectRaster(raster1,raster2)#将raster1的投影信息转换成与raster2一致
newproj <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +datum=WGS84" #自定义目标投影方式的WKT文本
RasterData <- projectRaster(raster1, crs=newproj)  #用crs参数指定投影信息

空间数据存储

栅格数据的导出

Terra包可以使用writeRaster函数输出栅格。

#raster是需要导出的栅格数据，filename是定义输出的路径和文件名
#format是输出栅格的格式，更多的格式可以参考https://www.rdocumentation.org/packages/raster/versions/3.5-29/topics/writeRaster
writeRaster(raster, filename=file.path(tmp, "test.tif"), format="GTiff", overwrite=TRUE)

矢量数据的导出

如果要将空间矢量数据以shp格式保存在本地，可以使用sf中的st_write()函数或它的别名函数wrtite_sf()。

st_write(area_wgs, "D:/Area_WGS84.shp")
write_sf(area_wgs, "D:/Area_WGS84.shp")

可视化制图(基于ggplot2)

在开始制图之前，我们需要先做一些准备
1. 制图数据的投影信息是否保持一致
2. 加载ggplot2、ggspatial和rasterVis库
3. 简单了解一些ggplot2中栅格数据的格式
下面是我们通过raster()函数读取进来的RasterLayer格式的栅格数据，数据结构大致如下图所示：

但这种RasterLayer格式却不能直接在ggplot2中进行绘制，需要转换成带X,Y坐标的DataFrame格式才能绘制,所以在可视化之前我们需要先把栅格数据从RasterLayer转换成df格式，然后才能进行绘图，详见下面的代码

#加载区域的线数据
boder <- read_sf("D:/Mydata/paper/data/粤港澳大湾地级市区线廓图.shp")
#对线数据进行坐标系转换
boder_wgs <- st_transform(boder,4326)
#将栅格数据转换成dataframe格式
df<- as.data.frame(as(mask_area,"Raster"),xy=T)
#自定义颜色
colors <-c("#33A02C","#B2DF8A","#FDBF6F","#1F78B4","#999999")
ggplot()+
#使用ggspatial包中的geom_tile()函数进行栅格数据可视化，用上图中的第三列进行fill填充
 geom_tile(data=df,aes(x=x,y=y,fill = Temp1990))+
 #比例尺设置
 annotation_scale(location = "bl") +
 #指北针设置
 annotation_north_arrow(location="tr",
                         style =north_arrow_nautical(
                           fill =c("grey40","white"),
                           line_col ="grey20"))+ 
  #将行政边界叠加在栅格数据上
  geom_sf(data=boder_wgs,aes(fill = NULL))+
  #栅格颜色设置，将NA值的颜色设置为transparent透明
  scale_fill_gradientn(colours=colors,na.value="transparent")+
  theme_bw()+
  #主题设置
  theme(panel.grid.major=element_blank(),
       panel.grid.minor=element_blank(),
       panel.background = element_blank(),
       legend.title = element_blank())

Raster包的函数介绍

Raster包是基于rgdal底层库编写的地理数据处理包，包括图像的创建，图像的读取，图像的运算，图像的可视化，图像的输出，非常实用

函数功能	函数名以及描述
增加一个栅格图层	addLayer(增加到的目标图层，增加图层1，增加图层2)
删除一个栅格图层	dropLayer(删除的目标图层，c(1，2，3)哪些层)
栅格栈	stack(图层1，图层2，图层3)
创建单个图层	raster(目标文件路径)
图层数	nlayers(目标)
找到目标位置周围的值	adjacent(目标图层, cells=哪些位置, directions=4周围邻域4 8 16, pairs=TRUE矩阵或者向量, target=NULL, sorted=F,include=F,id=F)
栅格合并 聚合函数	aggregate(x,fact=2,fun=mean)
一个限制对象（边界对象）	extent(xmin,xmax,ymin,ymax)
边界对齐，获得相同的原点和分辨率	alignExtent(e边界,r图层)
在不同图层之间切换，获得电影效果	amimate(X,pause=0.25,min,zlim,maxpixels=50000,n=10…)
产生n个随机数0-1	runif(n)
对图层中的NA进行插值	liner constant (X,filename=’’,method=“liner”,rule=1 or 2 描述如何处理第一个和最后一个)
计算 栅格面积 投影或者未投影	area()
栅格对象的算术方法	+ - * %% %/% 各图层直接运算。
转化为字符	as.character()
转化为数据框和矩阵	as.data.frame(X,row.names = 给行取名字,optional=T 以图层名字为列名,xy=T or F 坐标要不要)
建立一个图层列表对象	as.list(图层1，图层2)
转换逻辑矩阵 0为false 其他为true	as.logical(目标图层)
转换为矩阵或者向量，或者数组	as.matrix()
波段数，以及波段序数	r = raster(f,layer=2) nbands® bandnr®
栅格统计柱形图	barplot()
边缘检测	boundaries()
箱型图	boxplot()
创建栅格砖	brick() 可以由多种对象转化 extract(b,870)提取目标数据
对图层进行计算	calc(目标图层，fun 运算函数)
查找栅格值	cellFrom 很多种方法
通过边界查找	cellFromExtent(r,bb) bb边界
限制值的范围（任何高于上分位低于下分位都变成NA）	clamp(x,lower=,upper=)
清除内存中的栅格值	clearValues®
点击图层获得值	click®
增加等高线	contour(r,add=T)

原文链接：https://blog.csdn.net/qq_32832803/article/details/117257570 | |

本次的教程就差不多完结了，这里附上我看到的关于Raster包里一些函数的介绍，方便大家查看。还有其他关于技术交流的问题欢迎留言或者私信我~

参考文章
raster | R语言中的空间栅格对象及其基本处理方法（Ⅲ）：切片/掩膜、图层叠加
基于R语言进行栅格数据统计及Raster包简介
R语言Raster包和Terra包栅格读写、计算和一些使用经验分享
sf | 读取和保存空间矢量数据