学习日记 | R语言 | 多年NPP数据的读取和绘制均值时序变化图

最新推荐文章于 2024-04-27 17:33:21 发布

R ice

最新推荐文章于 2024-04-27 17:33:21 发布

阅读量1.2k

点赞数 4

文章标签：学习 r语言均值算法

本文链接：https://blog.csdn.net/Rigiser/article/details/128710150

版权

问题

从GEE上下载得到了多年的NPP（净生产力数据），想要对多年的栅格数据根据研究区域进行裁剪。这个工作可以在ArcGIS中直接通过【按掩膜提取】工具+【栅格统计】实现，实际操作中发现当栅格的波段数较多的时候，工作会变得繁琐，即使用上模型构建器，统计得到的均值仍需要逐一输入到表格中，再进行后续的绘图。而使用编程语言进行批量处理则效率更高。

数据来源

数据为GEE提供的MODIS NPP数据集，分辨率为500m。下载方式参考文章GEE下整洁地导出MODIS NPP和NDVI产品 - 知乎 (zhihu.com)

terra包

raster包和terra包都是R中处理栅格数据常用的包，综合对比和参考了其他文章后，发现terra的计算效率和方便性比raster包更高，但实际中还是需要两者都掌握，在不同情况下需要用到不同的包。

代码实现

首先批量读取文件夹中的栅格数据的路径名，list.files和list.dir都可以用来读取某个文件夹中的所有文件，pattern用于匹配名称具有特定文本的文件，$的作用是表示是以.tif结尾。

filename = list.files("F:/data",pattern=".tif$")
rasterpath = paste0('F:/data/',filename)

使用terra包中的rast函数读取栅格，vect函数用于读取矢量文件，读取完成后用plot绘制栅格数据查看。

library(terra)
npp=rast(rasterpath)
region = vect('F:/data/研究区域融合后.shp')
plot(region)
plot(npp)

在这里插入图片描述

如果使用raster包读取的话，代码如下，stack（）的作用是将多个向量整合成为栅格对象

library(raster)
npp = stack(rasterpath)

terra包和raster包读取栅格后会分别得到不同的栅格对象，terra包的栅格对象为spatraster，raster包的栅格对象为rasterlayer

使用terra包的mask进行掩膜，得到研究区域内的影像。

npp_clip = mask(npp,region)

使用global函数统计栅格数据的均值，得到一个列变量为’mean’的data.frame数据。

mean = global(npp_clip,'mean',na.rm=TRUE)

> mean2 #查看结果
          mean
a2003 1173.341
a2004 1141.611
a2005 1210.024
a2006 1189.618
a2007 1119.954
a2008 1152.369
a2009 1221.160
a2010 1221.160
a2011 1244.494
a2012 1333.663
a2013 1398.357
a2014 1340.123
a2015 1279.475
a2016 1350.586
a2017 1290.850
a2018 1425.772
a2019 1397.892
a2020 1364.782
a2021 1389.832

可以使用plot（）绘制时序变化点线图，$的作用是访问数据框中的’mean’列，type可以指定绘制的图的类型，这里选择了o（点线图）。

plot(seq(2003,2021),mean2$mean,type="o",pch=20,main="NPP均值变化",xlab="x轴",ylab="y轴")

在这里插入图片描述

用ggplot2也可以实现点线图的绘制。

![在这里插入图片描述](https://img-blog.csdnimg.cn/ccd41fad50ba4000a95256b95e48c0ac.png#pic_center)

library(ggplot2)
ggplot(data  = mean2,aes(x=seq(2003,2021),y=mean))
+geom_point(size=3)
+geom_line(size=0.8)
+xlab("年份")+ylab('NPP均值')