 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） ' e3 }$ g% d- c1 [+ R- c! D
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
: V4 x) r3 [4 M# ]9 K 
% X" C0 B* _0 i* j! O( s 本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
1 @: W7 ?5 E4 A 安装nppr包 . ~ l% F0 F' B' ?$ X+ Y
可在github（https://github.com/chaoxv/nppr）获取nppr包。
% b+ ^! S! o' Z+ v5 V+ m/ A #下载nppr包
/ _2 G/ k6 h+ Y/ G0 g+ S #install.packages(remotes)
$ o& M0 ]6 D; t8 L remotes::install_github(chaoxv/nppr) " |. ] ^3 h; j
#加载相关R包
G6 _6 r" g* l# w1 d* [ library(nppr)
4 U2 |$ x+ J' ?# M6 o9 Q library(RCurl) ! ^ x! x* Y1 l6 y0 h
library(XML) ; |* @# v9 r! _& m" C4 f' O
library(R.utils) $ D3 J0 Q1 J! h* W& t, [1 R
library(tidyverse)
+ B l( d5 R3 a7 Y4 ^, P: m) {: b( O library(lubridate) 3 f$ @* l" S% C2 I3 W
使用nppr包下载海洋遥感数据 ! \+ E. n: W! X/ }/ ]# V; \
nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 2 S0 e) l% R+ k% _; V
: \) {! I) Y- \8 R2 |) K N& U
以下以获取海洋NPP遥感数据为例作个演示。
3 b3 Z( a9 d* ^9 w0 b #创建工作路径
+ {# M$ T- {; Q/ c2 v) _8 i8 {! Q yourfolder <- F:/R/nppr/vgpm ' e: z5 e* }) H! }# V1 x: ^) K
dir.create(yourfolder)
, l% Z8 i2 h) G7 P, t# H #以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm
% v5 f: r% t `$ b2 l4 ] get_npp_vgpm( & A7 o/ j4 s F% x- {
file.path = yourfolder,
+ o7 X' A/ d/ T grid.size = low, #指定low或high可更改空间分辨率 $ K0 ?: p$ I0 Y1 j! N: R* q7 Y
time.span = monthly, #monthly代表月平均，dayly代表8天平均
- ?4 v5 _* h6 m- ?9 W6 a satellite = MODIS, #选择卫星
8 D3 \7 x1 S% Z( }; I mindate = 2016-01-15, #指定时间范围以下载遥感
* y3 [3 W3 R1 B- W% u* a4 @ maxdate = 2016-03-15
! d4 T! I9 y! o/ b5 z \7 k0 y9 b, F5 _ )
1 Z u" T2 i i1 I: L 
4 l$ N; \) T3 \4 ^* |
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
' u; C( y/ @0 G, J8 c% z- w6 k$ G9 C, B 使用nppr包进行遥感数据格式转换 / c6 \) d% a1 M8 t" {* r/ X6 Y
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 0 ^+ {! B! |4 r! [# _
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换 + C+ b" P6 n/ L+ ^4 U7 `) N
yourfile <- paste0(yourfolder, /201603.hdf)
0 p+ @" t! U9 l6 x vgpm <- read_hdf(file.path = yourfile)
; w- o/ v# p! k: z9 \ head(vgpm) 0 L, e* j; _4 |2 o! b
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
6 p9 l' H0 b# f2 |* L# L2 R 
4 A( z5 A- F+ t$ N v3 X4 ^. F 转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。 1 L4 P# y1 F: T" |" F! n) n
使用nppr包匹配目标经纬度的遥感数据
3 D% a; B6 Q5 ]0 V4 M 默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 5 G" j- }! B! t) P8 z+ N! U
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP 1 ^) D7 ^' i* b, l0 m, E6 {4 x
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
$ D% g& c! {5 E0 A4 }3 x2 }( m #或者同时指定多个数据，不再多说
7 T, h8 C- c# {5 B mydat <- data.frame( 7 [+ B2 x0 q" s1 i) d% l) q, r% [
lon = c(120, 112, 116),
8 y4 X# a& d3 b. p0 I lat = c(17, 15, 18), " x, y$ h, S/ |3 }" L$ `5 G( J
date = c(2016-03-04, 2016-03-07, 2016-02-04) # m( M! G( {. T' f Y
) & Q) q# b# W; w" m: s
match_df(mydat, file.path = yourfolder) 
1 N/ |# {' R/ u2 Z+ _$ M! p. S 绘制遥感地图 ; M' {( C5 o9 E% ?5 X
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
3 ]! q% H9 w7 _& j' F8 C #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 , s" }1 ]! w6 X3 B4 g! K
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP - J- E( B+ d% d( _' P) v
library(viridis) " [, x: k9 {5 y, ^" D
library(ggplot2)
% p0 l: C) F/ R8 M ggplot()+
7 t% u9 g% p: S1 n geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
c. `* M: _/ A" _7 A scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
$ g* B% D9 R! ?' w6 J' j labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
5 ?# j; Z M: [ 
$ C4 D5 Z8 {$ T0 ~1 {! y" F 根据时间和经纬度列表匹配遥感数据的批处理 1 P0 @# o9 _. I
实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） 2 F# g s; @& Z1 H3 c! ?- b
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 * G4 F, w8 ?: K
, B$ M( m `# D$ g2 T
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
! l' S7 g, D- r* L7 Y- S) G8 P ##如下以匹配SST数据为例做个演示 8 g/ A+ B# k9 N0 d0 | \8 ^
dat <- read.delim(data.txt) 7 S9 L3 `3 Y$ A2 N
Date <- unique(dat$Date) #获取日期 " @) [& _7 B! J. x9 w8 z K4 d
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 # Q$ B( M' p# J/ B
dir.create(yourfolder) " O) Q1 ^$ ~) f1 \5 {+ M! X% W& N# M
#通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
* _# M% W0 B- u2 ], j( l0 C4 D for (i in Date) { 6 s6 g7 K1 A1 X
yourfolder <- paste(getwd(), SST, i, sep = /) ( b! c; r$ y) u# d7 A
dir.create(yourfolder)
; j( K* h% V8 \+ ^% H# w get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
; N& X" w* I9 q yourfile <- dir(yourfolder)
/ ^) o2 H8 w" n8 f4 ~% q hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
& T$ u3 k }- z& R2 q6 ~/ n5 s write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) " M& P4 z( ?8 g2 n( }, j
}
. x3 L }: n2 {% u #再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均）
4 x. `7 n! c" o7 c3 b for (i in 1:nrow(dat)) {
! {5 T6 v* G8 X. J Date <- dat[i,Date] 9 K( L7 l( [) _6 J# c
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) ; B" w7 _& X- Z/ z) n* w$ v
hdf <- read.delim(yourfile) " p9 d: H3 ? R& u
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), ]
0 [6 D! V. I, P9 ` 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), ]
0 r8 B1 O2 G( t2 H0 r9 e }. |+ C dat[i,SST] <- mean(hdf$var) & }1 t+ Z9 d2 A! s# ?
}
5 t% ]2 m# @. r write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) . I! K2 m; U$ N! Y
" W+ }& @. i) w9 L9 W; ~2 s3 A
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
; C9 t4 j4 P; q) j, S% d8 C7 R% S/ C" c K1 U
4 N0 F$ G9 q+ b
3 @5 {* c! }& m- E3 I, z
1 H6 c, R6 V4 x+ G7 D0 O' z
| 
