R语言杂谈:给粉丝处理CMIP6的NC数据|提取东三省2015-2100平均温度降水
前几天有同学问我按照我的代码跑不出来,我跑的时候从没见过那个报错,我也很好奇为什么,然后就和他说等我最近忙完毕业论文再给他看看,然后最近正好提交系统外审,闲下来就跑了一下他的数据。解决完后发现有几点问题可以记录一下。
他的数据是CMIP6(做遥感的同学肯定不陌生),目的是提取东三省2015-2100的平均温度降水。

这个报错看起来是在导出栅格的时候,出现了一对多的问题。
后面我继续和这位同学沟通,发现他是一个大四的非农林专业的同学,对R语言以及遥感相当不熟练,无奈之下我只能问他可否把东三省的矢量和CMIP的nc数据发我一份我跑一下试试
(其实一开始我是想自己下,尽量避免拿别人的数据,但是CMIP下载实在是太慢了,东三省矢量我还得从全国shp里提取,也怕和他的研究区不一致,为了节约时间,只能让他发我了)

拿降水的nc数据举例

从这个基本信息中我们可以得知:这是一个拥有1032个图层(nlyr是n个layer是缩写)的栅格,储存着从2015年1月到2100年12月的逐月降水(Precipitation),注意它的单位是kg/m2/s。

关于坐标系的内容可以查阅这篇专栏。

然后在导入他的研究区矢量范围的时候发现居然没有坐标系信息


所以只能把shp定义一个坐标系,为了方便期间,直接用nc的坐标系
然后再看看nc栅格和shp的范围,探索一下


下一步就是按照矢量范围提取栅格均值
这一行的意思是,把pr这个nc数据先裁剪成东三省shp的外接矩形大小(用terra包的crop函数),再把这个外接矩形掩膜呈shp的形状(用terra包的mask函数),再提取这个shp形状的1032个图层的均值。最后会得到一个1032行的数据框,然后把它命名为df_pr。
(%>%是管道符号,表示前面的结果可以作为后面函数的输入)
至于为什么要先crop再mask可以参考下面的专栏。


通过head函数,我们可以看到基本上没有什么错误,可以进一步对数据框进行操作。
目前这个数据框一共有1032行,1列,行代表特定的年月,比如第一行pr_1代表2015年的1月,最后一行pr_1032代表2100年的12月,列表示东三省的降水均值,所以我们可以给她添加一列“年”以及一列“月”,来补充信息。

然后用openxlsx包的write.xlsx函数导出为xlsx

总结:
不要给我发任何关于你自己的涉密数据!!!尤其是你自己做实验得来的数据,但公开下载的数据可以发给我以便节省时间(比如MODIS,CMIP6这种谁都可以下载的)。因为我自己在学习初期请教别人的时候就发了好多我自己的关于我论文的数据,有的重要有的不重要,但是后来我小论文一直进展不顺利,迟迟没有发表,等我想投的时候发现已经有类似的了,作者和我当初请教的那个人是一个单位,甚至研究区和数据来源都一模一样,处理方法也和我的一致,很难不让人多想。由此导致的后果就是我不得不调整我的小论文结构,重新换题。
处理数据的过程中可能会有很多报错,可以先把报错内容搜索一下,探究一下是什么样的原因。