目录
- 0.问题导入
- 1.示例数据
- 2.导入NC数据并将其转化为data.frame
- 3.基于栅格的sc-PDSI的趋势项提取
- 4.EOF分析前数据结构调整
- 5.EOF分析
- 6.EOF模态解释方差分析（图1）
- 7.EOF时间系数可视化（图2）
- 8.EOF模态1-2可视化（图3）
- 9.总结
- 10.本篇所使用的软件包（没有的小伙伴需要用install.packages("")进行安装）
- 11.致谢
0. 问题导入
近些年EOF分析被广泛应用于研究气象要素的时空分布规律，那么这个操作如何实现呢？目前网上广为流传的是基于NLC与Matlab的，小编看了就头疼。。。于是乎有了今天这一篇，来基于R实现气象要素的EOF。
不出意外，本篇大概是网上仅有的R版EOF分析说明（中文版）
1. 示例数据
本篇所采用的NC示例数据与上一篇" R-ggplot2-说说绘图中颜色的这点事-应该是比较全的总结篇 "的示例数据一致，同样采用的是“全球sc-PDSI(干旱指数)1901-2018年的月尺度数据”。数据信息如下：
- 时间分辨率：月
- 空间分辨率：0.5°×0.5°
- 数据维度：360（行数）×720（列数）×1416（月份个数，也称数据深度），储存方式如图1。
下载链接如下：
点我下载NC文件
2. 导入NC数据并将其转化为data.frame
由于sc-PDSI 本身已经经过距平及标准化处理，本次EOF分析前不对数据进行以上处理。若是计划分析数据为降水，气温等，EOF分析前需要对数据进行距平及标准化的预处理。
setwd("/Users/jerseyshen/Documents/JianShu_Project/20191228/data")
ex = extent(100,120,40,48)
nc = stack("scpdsi_1901_2018.nc")
nc_crop = crop(nc,ex)
nc_df = as.data.frame(nc_crop,xy = T)
3. 基于栅格的sc-PDSI的趋势项提取
由于我们所采用的是月尺度的数据，其具有时间序列上的季节性。而EOF分析主要关注的是指标在时空尺度上演进的趋势，我们需要在分析前去掉指标（sc-PDSI）的内涵季节项，对趋势项进行提取。
nc_df_mat = nc_df[,-c(1,2)]
trend_decompose <- function(x){
ts_index = ts(x,start = c(1981,1),frequency = 12)
temp_trend = decompose(ts_index)$trend
na_index = which(is.na(temp_trend))
temp_trend = temp_trend[-na_index]
return(temp_trend)
}
nc_df_mat = apply(nc_df_mat,1,trend_decompose)
nc_df_mat = t(nc_df_mat)
4. EOF分析前数据结构调整
nc_df_mat = data.frame(long = nc_df[,1], lat = nc_df[,2],nc_df_mat)
na_index = c(1:6,635:640)
date = seq(as.Date('1901-01-01'),
as.Date('2018-12-01'),'1 month')
date = date[-na_index]
date = rep(date,each = nrow(nc_df_mat))
nc_df_melt = melt(nc_df_mat,c('long','lat'))
nc_df_melt$date = date
5. EOF分析
source('~/Documents/JianShu_Project/20191228/data/EOF_modified.R')
eof = EOF_modified(value ~ lat + long | date, 1:10,data = nc_df_melt,
B = 1000, probs = c(low = 0.1, hig = 0.9))
6. EOF模态解释方差分析（图1）
根据图1可见，模态1-2方差解释率分别为41%与20%，一般认为一般认为方差解释率小于 10%的气候要素场的气候物理意义较小。因此，本示例选择前两个模态进行分析。
summary(eof)
Importance of components:
Component Explained variance Cumulative variance
1 41% 41%
2 20% 61%
3 10% 71%
4 6% 76%
5 4% 81%
6 3% 83%
7 2% 86%
8 2% 87%
9 1% 89%
10 1% 90%
p1 = screeplot(eof)+theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
axis.title = element_text(face = 'bold',color = 'black',size = 14,hjust = 0.5),
legend.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.height=unit(0.1,"inches"),
legend.key.width=unit(0.5,"inches")
)+
xlab('PC')+
ylab('R-squared')
png('plot1.png',
height = 20,
width = 20,
units = 'cm',
res = 800)
print(p1)
dev.off()
图1 EOF各模态解释方差分析
7. EOF时间系数可视化（图2）
aao1 = cut(eof,1)
aao2 = cut(eof,2)
pl_line_df = rbind(aao1$right,aao2$right)
p2 = ggplot()+
geom_line(data = pl_line_df,aes(x = date, y = value,color = PC))+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
axis.title = element_text(face = 'bold',color = 'black',size = 14,hjust = 0.5),
legend.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.height=unit(0.1,"inches"),
legend.key.width=unit(0.5,"inches")
)+
xlab('PC')+
ylab('Indices')
png('plot2.png',
height = 10,
width = 20,
units = 'cm',
res = 800)
print(p2)
dev.off()
图2 EOF时间系数可视化
8. EOF模态1-2可视化（图3）
pl_spatial_df = rbind(aao1$left,aao2$left)
pl_spatial_df$cuts = cut(pl_spatial_df$value,
breaks = seq(-0.09,0.08,0.02))
p3 = ggplot()+
geom_tile(data = pl_spatial_df,aes(x = long,y = lat, fill = cuts))+
scale_fill_brewer(palette = 'Spectral')+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
axis.title = element_text(face = 'bold',color = 'black',size = 14,hjust = 0.5),
legend.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.height=unit(0.1,"inches"),
legend.key.width=unit(0.5,"inches")
)+
facet_wrap(~ PC, ncol = 2)+
xlab('Longitude')+
ylab('Latitude')
png('plot3.png',
height = 10,
width = 20,
units = 'cm',
res = 800)
print(p3)
dev.off()
图3 EOF模态1-2可视化
9. 总结
本篇主要解决了以下三个问题：
- 如何在R中完成时空EOF分析
- 如何从EOF分析中选择模态？
- 如何在将EOF分析结果可视化？
对于本篇结果的讨论，小编在此处挖个坑，等待各路大神来在评论区内讨论～
10. 本篇所使用的软件包（没有的小伙伴需要用install.packages("")进行安装）
特别说明：metR中的EOF分析函数有BUG，小编已Debug～如有需要烦请转发本篇博文并集够10个赞，然后联系小编获取哈～
library(metR)
library(reshape2)
library(raster)
library(ncdf4)
library(viridis)
library(ggplot2)
11. 致谢
感谢大家的持续关注，小编会继续努力，持续更新下去的！
大家如果觉得有用，还麻烦大家转发点赞加关注哈，也可以扩散到朋友圈，多谢大家啦～
大家如果在使用本文代码的过程有遇到问题的，可以留言评论，也可以私信我哈~~
小编联系方式
网友评论