R-如何在NC时间堆叠数据集中根据坐标提取数据并展示其趋势

0.问题导入
1.示例数据
2.导入NC数据与兴趣点坐标
3.根据兴趣点左边范围裁剪shp文件
4.将兴趣点转化为SpatialPoints
5.根据兴趣点坐标信息提取对应位置sc-PDSI指标时间序列
6.sc-PDSI序列矩阵NA值检查
7.计算兴趣点对应序列的趋势值
8.空间可视化
9.空间可视化优化
10.本文所用到的软件包
11.本文总结
12.致谢

0. 问题导入

当我们拿到一个高维的NC数据的时候，可能有时候不需要对所有栅格值进行计算，只是单纯地想根据兴趣点提取对应位置的指标时间序列，那么我们应该怎样处理呢？本篇文章给出解决方案。

同之前一样，本文所用软件包附于文末

1. 示例数据

本篇所采用的NC示例数据与上一篇" R-ggplot2-说说绘图中颜色的这点事-应该是比较全的总结篇 "的示例数据一致，同样采用的是“全球sc-PDSI(干旱指数)1901-2018年的月尺度数据”。数据信息如下：

时间分辨率：月
空间分辨率：0.5°×0.5°
数据维度：360（行数）×720（列数）×1416（月份个数，也称数据深度），储存方式如图1。

下载链接如下：
点我下载NC文件

图1 NC文件数据储存格式－吃货解说版

本篇采用示例坐标数据为中国区域内随机生成的坐标,下载方式如下:
点我下载兴趣点数据

本篇采用的世界大洲级地图(shp文件)下载方式如下:
点我下载world.shp 文件

2. 导入NC数据与兴趣点坐标

input_nc = 'L:\\JianShu\\2019-12-07\\data\\scpdsi_1901_2018.nc'
nc_stack = stack(input_nc)

input_loc = 'L:\\JianShu\\2019-12-08\\data\\loc.csv'
loc = read.csv(input_loc,header = T)
loc = loc[,-1]

3. 根据兴趣点左边范围裁剪shp文件

ex_crop = extent(min(loc[,1])-2,max(loc[,1])+2,min(loc[,2])-2,max(loc[,2])+2)
world = shapefile('L:\\JianShu\\2019-12-08\\data\\shp\\world_continent2.shp')
world_crop = crop(world,ex_crop)

4. 将兴趣点转化为SpatialPoints

sp = SpatialPoints(loc)

5. 根据兴趣点坐标信息提取对应位置sc-PDSI指标时间序列

df = extract(nc_stack,sp)
df = t(df)
print(ncol(df)) #列数对应为点的个数

6. sc-PDSI序列矩阵NA值检查

col_na = which(is.na(apply(df,2,mean)))
if(length(col_na) == 0){
    df_na_re = df
    loc_na_re = loc
}else{
    df_na_re = df[,-col_na]
    loc_na_re  = loc[,-col_na]
}

7. 计算兴趣点对应序列的趋势值

趋势值的计算我们采用的是改进后的Mann-Kendall 趋势分析方式, 详细介绍可参照R-Modified Mann-Kendall 趋势分析（改进后MK趋势分析）这篇博文, 里边文末有附本文所用的mmkTrend函数, 记住一定将其保存为R文件,然后在Rstudio中source一下,才能用哈~就像这样:

 source('L:/JianShu/2019-12-08/mmkTrend.R')

mmk_cal<-function(x){
  temp_mk = mmkTrend(x)$Zc
  return(temp_mk)
}
mmk_box = apply(df_na_re,2,mmk_cal)

8. 空间可视化

df = data.frame(
long = loc_na_re[,1],
lat = loc_na_re[,2],
MMK = as.numeric(mmk_box)
)
df_m = melt(df,c('long','lat'))

df_m$cuts = cut(df_m$value,
breaks = seq(round((min(df_m$value)-1)),round((max(df_m$value)+1)),1))

df_m$size = df_m$cuts

world_crop_df = fortify(world_crop)

p = ggplot()+
  geom_polygon(data = world_crop_df,aes(x = long,y = lat, group = group),
               fill = 'transparent',color = 'black',size = 0.5)+
  geom_hline(aes(yintercept = 50),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_hline(aes(yintercept = 40),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_hline(aes(yintercept = 30),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_hline(aes(yintercept = 20),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_vline(aes(xintercept = 80),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_vline(aes(xintercept = 100),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  geom_vline(aes(xintercept = 120),linetype = 'dashed',alpha = 0.5,lwd = 0.5,color = 'black')+
  
  geom_point(data = df_m,aes(x = long,y = lat, color = cuts,size= size))+
  theme(panel.background = element_rect(fill = 'transparent',color = 'black'),
        axis.text = element_text(face='bold',colour='black',size=fontsize,hjust=.5),
        axis.title = element_text(face='bold',colour='black',size=fontsize,hjust=.5),
        legend.position=c('bottom'),
        legend.direction = c('horizontal'))+
  scale_color_brewer(palette = 'Spectral')+
  scale_size_manual(values = 1:10)+
  coord_fixed(1.3)+
  
  guides(fill=guide_legend(nrow=2))+
  xlab('Longitude')+
  ylab('Latitude')

图2 兴趣点1901-2018年sc-PDSI MMK趋势空间分布图

但是，大家有没有觉得很别扭

嗯嗯，是的！！！又是一个逼死强迫症患者的技术BUG，点大小的分级与颜色的分级没有重合！！！是不是很难受！

来, 下边我们一起来解决这个BUG!

9. 空间可视化优化

首先, 问题诊断. 出现以上问题的原因呢是在这一句:

 geom_point(
data = df_m,
aes(x = long,y = lat, 
color = cuts,
size= size))+

很多时候,我们为了绘图理解的方便,会将color和size的分组在绘图的data.frame中分开,尽管他们是一样的。

head(df_m)
    long   lat variable      value    cuts    size
1 122.52 52.97      MMK  1.8868017   (1,2]   (1,2]
2 124.48 47.80      MMK  0.9688879   (0,1]   (0,1]
3 125.25 45.70      MMK -0.4857201  (-1,0]  (-1,0]
4  84.67 44.43      MMK  1.1515050   (1,2]   (1,2]
5  75.25 39.72      MMK  2.2346846   (2,3]   (2,3]
6 101.07 41.95      MMK -2.3940473 (-3,-2] (-3,-2]

但这正是导致问题的症结所在, 为了将两组图例合并在一起,我们应该将出问题的这一句绘图函数更换成下式:

 geom_point(
data = df_m,
aes(x = long,y = lat, 
color = cuts,
size= cuts))+

然后,绘图结果就如图3所示, 将size与color的图例合并为一个,这样看起来就舒服多了~~~

图3 优化后的图件

10. 本文所用到的软件包(没有的记得用install.packages('package name')进行安装)

library(raster)
library(reshape2)
library(sp)
library(RColorBrewer)
library(ggplot2)

11. 总结

回顾一下, 本篇文章主要解决如下几个问题：

如何从NC时间堆叠数据集中根据坐标提取数据?
如何计算各个兴趣点的MMK趋势值？
如何根据各点MMK趋势值的大小绘制空间图并体现其差异性？
当点的大小与颜色的图例分开时，如何将其合并在一起？

12. 致谢

大家如果觉得有用，还麻烦大家关注点赞，也可以扩散到朋友圈，帮助到绘图中同样陷入颜色选择困难症的TA

大家如果在使用本文代码的过程有遇到问题的，可以留言评论，也可以私信我哈~~

我的联系方式