目錄
- 0.問(wèn)題導(dǎo)入
- 1.示例數(shù)據(jù)
- 2.導(dǎo)入NC數(shù)據(jù)與興趣點(diǎn)坐標(biāo)
- 3.根據(jù)興趣點(diǎn)左邊范圍裁剪shp文件
- 4.將興趣點(diǎn)轉(zhuǎn)化為SpatialPoints
- 5.根據(jù)興趣點(diǎn)坐標(biāo)信息提取對(duì)應(yīng)位置sc-PDSI指標(biāo)時(shí)間序列
- 6.sc-PDSI序列矩陣NA值檢查
- 7.計(jì)算興趣點(diǎn)對(duì)應(yīng)序列的趨勢(shì)值
- 8.空間可視化
- 9.空間可視化優(yōu)化
- 10.本文所用到的軟件包
- 11.本文總結(jié)
- 12.致謝
0. 問(wèn)題導(dǎo)入
當(dāng)我們拿到一個(gè)高維的NC數(shù)據(jù)的時(shí)候谣妻，可能有時(shí)候不需要對(duì)所有柵格值進(jìn)行計(jì)算，只是單純地想根據(jù)興趣點(diǎn)提取對(duì)應(yīng)位置的指標(biāo)時(shí)間序列，那么我們應(yīng)該怎樣處理呢僚匆？本篇文章給出解決方案。
同之前一樣韭山，本文所用軟件包附于文末
1. 示例數(shù)據(jù)
本篇所采用的NC示例數(shù)據(jù)與上一篇" R-ggplot2-說(shuō)說(shuō)繪圖中顏色的這點(diǎn)事-應(yīng)該是比較全的總結(jié)篇 "的示例數(shù)據(jù)一致箫章，同樣采用的是“全球sc-PDSI(干旱指數(shù))1901-2018年的月尺度數(shù)據(jù)”。數(shù)據(jù)信息如下：
- 時(shí)間分辨率：月
- 空間分辨率：0.5°×0.5°
- 數(shù)據(jù)維度：360（行數(shù)）×720（列數(shù)）×1416（月份個(gè)數(shù)系谐，也稱(chēng)數(shù)據(jù)深度），儲(chǔ)存方式如圖1讨跟。
下載鏈接如下：
點(diǎn)我下載NC文件
本篇采用示例坐標(biāo)數(shù)據(jù)為中國(guó)區(qū)域內(nèi)隨機(jī)生成的坐標(biāo),下載方式如下:
點(diǎn)我下載興趣點(diǎn)數(shù)據(jù)
本篇采用的世界大洲級(jí)地圖(shp文件)下載方式如下:
點(diǎn)我下載world.shp 文件
2. 導(dǎo)入NC數(shù)據(jù)與興趣點(diǎn)坐標(biāo)
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. 根據(jù)興趣點(diǎn)左邊范圍裁剪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. 將興趣點(diǎn)轉(zhuǎn)化為SpatialPoints
sp = SpatialPoints(loc)
5. 根據(jù)興趣點(diǎn)坐標(biāo)信息提取對(duì)應(yīng)位置sc-PDSI指標(biāo)時(shí)間序列
df = extract(nc_stack,sp)
df = t(df)
print(ncol(df)) #列數(shù)對(duì)應(yīng)為點(diǎn)的個(gè)數(shù)
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. 計(jì)算興趣點(diǎn)對(duì)應(yīng)序列的趨勢(shì)值
趨勢(shì)值的計(jì)算我們采用的是改進(jìn)后的Mann-Kendall 趨勢(shì)分析方式, 詳細(xì)介紹可參照R-Modified Mann-Kendall 趨勢(shì)分析（改進(jìn)后MK趨勢(shì)分析）這篇博文, 里邊文末有附本文所用的mmkTrend函數(shù), 記住一定將其保存為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')
但是纪他，大家有沒(méi)有覺(jué)得很別扭
嗯嗯， 是的Ａ澜场２杼弧！又是一個(gè)逼死強(qiáng)迫癥患者的技術(shù)BUG凉馆，點(diǎn)大小的分級(jí)與顏色的分級(jí)沒(méi)有重合Ｐ皆ⅰＭ鲎省！是不是很難受向叉！
來(lái), 下邊我們一起來(lái)解決這個(gè)BUG!
9. 空間可視化優(yōu)化
首先, 問(wèn)題診斷. 出現(xiàn)以上問(wèn)題的原因呢是在這一句:
geom_point(
data = df_m,
aes(x = long,y = lat,
color = cuts,
size= size))+
很多時(shí)候,我們?yōu)榱死L圖理解的方便,會(huì)將color和size的分組在繪圖的data.frame中分開(kāi),盡管他們是一樣的锥腻。
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]
但這正是導(dǎo)致問(wèn)題的癥結(jié)所在, 為了將兩組圖例合并在一起,我們應(yīng)該將出問(wèn)題的這一句繪圖函數(shù)更換成下式:
geom_point(
data = df_m,
aes(x = long,y = lat,
color = cuts,
size= cuts))+
然后,繪圖結(jié)果就如圖3所示, 將size與color的圖例合并為一個(gè),這樣看起來(lái)就舒服多了~~~
10. 本文所用到的軟件包(沒(méi)有的記得用install.packages('package name')進(jìn)行安裝)
library(raster)
library(reshape2)
library(sp)
library(RColorBrewer)
library(ggplot2)
11. 總結(jié)
回顧一下, 本篇文章主要解決如下幾個(gè)問(wèn)題：
- 如何從NC時(shí)間堆疊數(shù)據(jù)集中根據(jù)坐標(biāo)提取數(shù)據(jù)?
- 如何計(jì)算各個(gè)興趣點(diǎn)的MMK趨勢(shì)值？
- 如何根據(jù)各點(diǎn)MMK趨勢(shì)值的大小繪制空間圖并體現(xiàn)其差異性母谎？
- 當(dāng)點(diǎn)的大小與顏色的圖例分開(kāi)時(shí)瘦黑，如何將其合并在一起？
12. 致謝
大家如果覺(jué)得有用奇唤，還麻煩大家關(guān)注點(diǎn)贊幸斥，也可以擴(kuò)散到朋友圈，幫助到繪圖中同樣陷入顏色選擇困難癥的TA
大家如果在使用本文代碼的過(guò)程有遇到問(wèn)題的咬扇，可以留言評(píng)論甲葬，也可以私信我哈~~
