使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP)
1 J% k/ [# X! w Ocean Productivity(http://sites.science.oregonstate.edu/ocean.productivity/index.php)是一个众所周知的海洋生产力数据库,我们经常从中下载相关的遥感数据来用于分析。 9 a7 e3 h s. a3 }
 - E+ X0 E' q3 m+ R& o2 X
本篇介绍师兄的一个R包,nppr包(https://github.com/chaoxv/nppr)。该包提供了便捷的函数,可以用来下载和处理Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。
* O0 V( R; U6 v7 d$ Q- k2 Y 安装nppr包
; l* O4 ~8 w, a! ~) I1 u! I9 O$ \ 可在github(https://github.com/chaoxv/nppr)获取nppr包。
- f5 G6 i/ ~/ p0 a: P( D #下载nppr包
5 G4 a* W2 w% H1 O6 b! n9 b #install.packages(remotes)
+ G: P- W+ Q2 u# r4 |, `/ {$ k remotes::install_github(chaoxv/nppr)
7 v" ?5 f4 g6 C+ g; k #加载相关R包
$ _2 x2 x/ m+ N( W+ I$ _% k library(nppr) ?- z; I p" `! M
library(RCurl) ! a0 f! V' U, G' U: P6 o( A
library(XML) , _" n3 ~2 S% v3 O- l9 R$ C
library(R.utils) W8 s3 Y8 U8 v- s' s) W
library(tidyverse)
: Y2 `4 f( T7 g# Q0 R/ h" j library(lubridate) 1 J+ u1 ?; q' q) Y. t- V9 [
使用nppr包下载海洋遥感数据
! P3 N# P$ u- w: F nppr包内常用函数如下所示,get_*等函数可分别用于下载Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。
2 t, t4 R. G' {- {- t6 @  - l. {: s) A. a2 S
以下以获取海洋NPP遥感数据为例作个演示。 a, x; ^$ g+ n2 W5 o& u7 v9 ~ c
#创建工作路径 9 }# }* ~& L) f
yourfolder <- F:/R/nppr/vgpm
8 g6 V5 w ?8 D5 } dir.create(yourfolder)
! S6 ^4 o2 P5 |. U k7 h #以VGPM(NPP的一种遥感算法)为例做个演示,详情?get_npp_vgpm
5 Z* I% }$ I: T l4 ?5 A get_npp_vgpm( 8 p G- k! d) C* z2 N! S
file.path = yourfolder, ) U$ D) S8 w: \$ F3 X) V, x6 j
grid.size = low, #指定low或high可更改空间分辨率
/ E9 I; {; n: |7 i8 J+ P time.span = monthly, #monthly代表月平均,dayly代表8天平均
- i! K4 C2 a" d I7 A7 _ satellite = MODIS, #选择卫星 ) N# a; J3 p0 O+ C1 M( z
mindate = 2016-01-15, #指定时间范围以下载遥感
v3 h& o. }0 H# b, E( o) O) A6 h maxdate = 2016-03-15
# l% T3 c3 h- z6 `/ R- l3 \ ) - v4 l' p7 V! p6 ^) o
 0 O$ p1 p" h; U- N
在这个示例中,我们使用nppr包下载了来自Ocean Productivity(http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php)的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
+ W: ]' }7 |1 j1 \ D x8 a 使用nppr包进行遥感数据格式转换 . A5 k2 H# a7 Z' a) l& P1 X
如上所述,下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法,可将hdf格式转为常见的数据框格式,便于我们后续操作。 & d. X6 @1 q1 A! I& \! T
#将hdf文件转为常见的数据框格式,例如将上述下载的2016年3月的月平均NPP数据做个转换
( ]0 F6 H# e+ c9 Q- u yourfile <- paste0(yourfolder, /201603.hdf)
# y" |& [8 C6 t0 D vgpm <- read_hdf(file.path = yourfile)
2 P2 e. _; N& e) E* s# E+ E+ c head(vgpm)
9 N1 }' U1 h5 e write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
V' h$ w( W; A( k7 P  " L# M- k9 j& l6 @4 r s
转换后的数据框包括三列,分别是经度(lon)、纬度(lat)以及当前日期内该经纬度海区的NPP(var,单位mg C m-2 d-1)。 . D: c( K7 H$ e: i3 X7 E" ?0 x
使用nppr包匹配目标经纬度的遥感数据 & X# c( X3 [' b
默认情况下,下载的遥感数据是全球海洋的。nppr包同样提供了相关函数,便于我们从中提取特定区域的遥感数据,如下所示。 9 \ }3 k/ E* f; L
#获取指定日期和经纬度的遥感,例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
% o% H( H3 q+ j z' @ match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
6 ]; H7 m2 L9 r; w% ?: J# U W( f #或者同时指定多个数据,不再多说
3 w/ f. U) y# u) b& u& Q4 M$ D- Z; C3 w mydat <- data.frame( : i8 ^6 g- s5 R( u9 K4 N+ H
lon = c(120, 112, 116),
8 C( \; p" Y p: F( m2 Y5 R4 h lat = c(17, 15, 18),
( U O' m! r" c7 N4 d9 F4 ]+ B date = c(2016-03-04, 2016-03-07, 2016-02-04) 1 U# E, Q! _6 t" L: U% i
)
B4 |0 a0 C0 k M# S3 L match_df(mydat, file.path = yourfolder) ( d- [" t/ o+ H2 r8 d% v
绘制遥感地图
: e+ C! y( B( C nppr包的函数geom_oce()可以用来绘制地图,例如我们作图展示来自遥感反演的NPP分布。
0 Z( @! k! {& D #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 ; n8 o) s1 s% l* R
#我们仍以该数据为例作图,展示中国南海2016年3月的月平均NPP / ~# w% o$ U, }! F
library(viridis) & Z3 U, e( P7 d: o6 S R' e+ M" Y2 I
library(ggplot2)
7 ~7 D% e' V b) E) v: I, s, j ggplot()+ 1 H4 M; h8 b+ ~" c
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ , _& f9 X! _' Q: `( i( y3 `# y- j8 e
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ 4 |: v1 G, Y' h& o9 u# b5 Q b
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
. }2 ^$ B! y/ A Z 
- g5 C7 f. R- u8 r& } 根据时间和经纬度列表匹配遥感数据的批处理
- Y; _" c( F% X. \0 e" X4 m- } 实际情况中,经常需要对来自不同时间不同经纬度的大量站位匹配遥感,以下提供了一个批处理(不过这是自己先前瞎写的,然后一直偷懒一直用,俺也不知道写的对不对......写在这里仅为方便自己复制粘贴,大家慎用......)
. r3 e( x* j0 a1 @! u$ q3 M 将待匹配的站位的经纬度、日期信息整合在一个文档中,如下所示的这样(本示例命名为“data.txt”)。
) y% F) F3 G6 I$ x$ |/ F. t 
8 S3 X- V. ^8 }& M2 \ 随后在R中读取该文件,设置一个循环,依次读取日期信息以下载当前日期的遥感(如月平均或8天平均的SST、PAR、Chla或NPP等)。并再根据各站位的经纬度,从中匹配该站位附近的数据(比方说以0.1°为网格进行匹配,并将网格内的数据平均)。
6 j& J$ y# r+ Y5 c) K1 f ##如下以匹配SST数据为例做个演示
7 ~! M% r% W2 B dat <- read.delim(data.txt)
8 I- L( D( f4 t: j Date <- unique(dat$Date) #获取日期
; ]# R, @! P8 O yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 2 i/ q3 z* u8 q
dir.create(yourfolder) & _/ j* ], [- L: H/ t
#通过循环依次获取各日期下的遥感(本示例以下载8天平均SST为例)
2 D( i: O! x) t. K' n, d for (i in Date) { 9 v8 {/ N! ]4 B6 s1 B& y( H' B
yourfolder <- paste(getwd(), SST, i, sep = /) 3 \ T: S4 @9 w1 a: C: k: d
dir.create(yourfolder) 0 ]! B+ e1 A+ X! J5 ^( S6 P
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) % O; h) g2 W3 I- r R$ q# ]7 R# w
yourfile <- dir(yourfolder) : M: f( E# F% ? x- C4 g, b# H
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) 1 g& L5 n* W6 ^8 `1 r
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
! t( p0 A' H& R: \0 r8 i }
" a& E B; W! I8 ?# b #再根据列表中各站位的经纬度匹配当前日期的遥感(本示例计算0.1°网格内的平均) : v) D* W. P) V" X. ]5 Z, @
for (i in 1:nrow(dat)) {
% T5 Y k) Z' | Date <- dat[i,Date]
# Z$ e) r, m6 a$ } yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) ; T4 c4 Q f5 w. A
hdf <- read.delim(yourfile) . |8 D; c* z; D2 T( r
hdf <- hdf[which(round(hdf$lon, 2) < round(dat[i,Longitude], 2)+0.1 & round(hdf$lat, 2) < round(dat[i,Latitude], 2)+0.1), ]
9 [( `5 \9 V( h hdf <- hdf[which(round(hdf$lon, 2) > round(dat[i,Longitude], 2)-0.1 & round(hdf$lat, 2) > round(dat[i,Latitude], 2)-0.1), ]
8 z: ~: ~' \6 ~8 x; t. g( { dat[i,SST] <- mean(hdf$var)
/ y# \5 U+ o' _1 ?; B/ H, _. K }
2 G! N4 y# ^+ N1 Y J4 c3 x# A write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
, x. C( O$ h2 Q 
2 g' l% |$ ]$ }) p9 Y& I# D 输出列表的最后一列添加了匹配的遥感数据(本示例匹配了SST)。
& T# U& V4 f! P1 L2 @/ P% }- S
8 Z. e1 f* q( I5 r+ O" R; r. G# V6 k
$ Y: E' Y3 i" u. {
% }# e0 e# I- A8 @* ^( \3 e
|