 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） $ D. V9 r$ t# s2 p& L# a
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
! C6 Y: p2 {* B. V 
% h- z" Q4 ]4 F# B5 k) W
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。  5 d M6 L0 n& X( G
安装nppr包 - }* f% y5 P! P( P7 n% }) _; B
可在github（https://github.com/chaoxv/nppr）获取nppr包。
4 B( [9 f" s' x #下载nppr包
5 [& S" k4 ?7 g5 F #install.packages(remotes) 1 a0 F/ }' W$ q# o( H8 g
remotes::install_github(chaoxv/nppr) , \# j* [# u5 v' S* G) c( {
#加载相关R包
( x% M( T( }0 I4 _ B3 ~8 N1 H library(nppr) ! v& s+ X e. K, {: y1 L* r
library(RCurl)
# m7 l/ h* p, v library(XML)
) I5 ]* B( `, v; X' I: o2 w0 y library(R.utils)
% u7 s' n. \2 k0 ] library(tidyverse)
5 x: N( z$ {4 y6 F library(lubridate)
% V9 ~( u7 }) ^1 N8 C6 Q) { 使用nppr包下载海洋遥感数据
' \+ u9 s) x7 p! ?- W- j# [ nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 " j: c1 H2 J) E: [" W
! q2 w6 e( L7 C+ b* K3 o
以下以获取海洋NPP遥感数据为例作个演示。
, @' V( Q5 g( S& t/ R; r( e* u4 n. F* M #创建工作路径
7 \2 W2 \7 j4 a- [! B yourfolder <- F:/R/nppr/vgpm # D5 P4 h' C9 y' w: X# U# p' x
dir.create(yourfolder) % G% l2 L5 u4 {/ [) ]$ ?
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm
% d2 ], j9 T8 P! {6 p: m' J get_npp_vgpm( 4 `4 |" p+ x& i% d; e8 X9 w: I" w
file.path = yourfolder,
0 f/ r V& N; ~& u/ B grid.size = low, #指定low或high可更改空间分辨率
; T6 m# s4 Q8 r* I0 ^2 B8 c time.span = monthly, #monthly代表月平均，dayly代表8天平均
' O# t, w$ |9 K% [) q O satellite = MODIS, #选择卫星 9 L2 t- S9 Y4 ?1 @3 y
mindate = 2016-01-15, #指定时间范围以下载遥感
5 w3 N3 ~1 {, Y maxdate = 2016-03-15 ) \: B [7 z' v/ b* l: p
) 7 T$ L, m( v1 t- x5 }" `* |+ G
- \% O( J$ i: w! S' k
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 5 l) O) k" V9 F z1 G3 [; P0 f0 u
使用nppr包进行遥感数据格式转换
, Z# f2 n2 I Z; o# a! k$ P9 o* g 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 : P% s W( ]( J& m `, J' e
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
8 R5 c* k7 u* Z yourfile <- paste0(yourfolder, /201603.hdf) # z2 Y: O% \& y$ B- ^1 Y
vgpm <- read_hdf(file.path = yourfile) # w% e* t0 L7 \; ?) x: H
head(vgpm)
8 {! C/ E2 i2 { Z write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) & b6 D+ q( W- {) w. o7 j% s
* y7 j; l3 a1 ?2 D+ Y0 r& g 转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
7 Z- l! v1 k! ?8 x6 c 使用nppr包匹配目标经纬度的遥感数据 : l2 r4 K4 N T+ y( z( \+ R
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。
5 h/ N! q+ Q: ?, q/ O #获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP ' e7 C4 I1 R o6 e, S
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) $ L) m1 u5 }* A) {. ?) T1 ]
#或者同时指定多个数据，不再多说 8 y/ R9 s, [, N3 l
mydat <- data.frame( 1 g" E/ `; ~/ O3 w9 s2 j
lon = c(120, 112, 116), - I- |2 G/ B8 Q7 `# |) @0 R
lat = c(17, 15, 18),
4 M; Y w5 W& ?. I/ M date = c(2016-03-04, 2016-03-07, 2016-02-04) * _% n Q( c n4 F! e$ M2 g
) 9 x9 W8 r3 Q* Q2 J) ~
match_df(mydat, file.path = yourfolder)  . L; E& @; L0 D( J# B% n
绘制遥感地图
1 B0 v* J- O, J2 \ nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
( V( M" D) M* k2 {& D #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
3 \7 o4 v& d5 T6 | #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP 1 l; O7 B: R" l
library(viridis) ' O- ?# K: q2 A
library(ggplot2) 3 G! O+ _1 c: G: @* e& Z0 Q2 ]" A# K
ggplot()+
" q3 \- x, ~. y. X$ d: I geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
! e$ e4 i7 ` h scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
, U; f9 r W9 x. k+ ? labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) 8 H+ R1 w; V0 l+ j
. m; e1 i' P6 s. U: g6 y( ^
根据时间和经纬度列表匹配遥感数据的批处理
- g/ o L' I0 P* T0 d: K 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
- \8 b- R _" D: }: D& [ 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 1 a' f5 M, A+ x8 b" j, J
) O- `9 j3 @2 B. d 随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。 1 ?0 [+ D. h" h! w) G4 o; G0 ^
##如下以匹配SST数据为例做个演示 : ^" y& b8 L' s. Q1 d, t- C
dat <- read.delim(data.txt)
" D' y3 b$ n, ~4 w) S6 O! G4 ^8 C# ?- q Date <- unique(dat$Date) #获取日期
- @3 H. [8 q6 \9 g1 l+ ], Q& i. A yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 ! { G m& n1 v Z( m
dir.create(yourfolder)
' F' R4 y! u0 l( m$ O6 J4 k [7 o #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
4 r( b( Q2 o) {* x4 O' q) `3 R for (i in Date) {
4 b' v5 N+ y0 [. D3 P! \- N yourfolder <- paste(getwd(), SST, i, sep = /) . ~$ H _! e9 p0 m# p1 M+ b
dir.create(yourfolder)
) n* B4 X) J5 l/ Y, u4 ^ get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
2 o. ]8 c$ `6 E" K8 G7 G yourfile <- dir(yourfolder)
# I" T9 Y) E- n/ l% ~ hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
" z7 @% R8 V7 H! Y: D1 V write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
) B. t4 O6 T9 F' V } * D$ c+ s4 W+ h# H; p$ B
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均）
( x- @' l, ?3 D) ~% a, Y( p for (i in 1:nrow(dat)) {
8 h( O+ K/ ~6 S: y Date <- dat[i,Date]
6 l; }9 Q; s% V4 U1 M yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) : C; V& k+ z1 l2 C
hdf <- read.delim(yourfile) " n" r! m# x0 p/ [
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), ]
& Y& l. V$ E# S% K& 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), ] # m" y9 C! l; F5 g
dat[i,SST] <- mean(hdf$var)
2 g! e+ U0 a5 N) O% I: M% R }
/ s8 F {+ \% t- c write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) 4 v! \& `1 r" S$ V
$ I/ H& d) ~% u! A
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
) a" o! I3 K1 f3 J3 R" U) D! r
6 O8 K, {4 R$ g+ ^. q$ ^3 M$ v4 x5 t! h- d' q
4 K) M( h( z5 \# p
8 F1 M, N5 w3 Z) T/ j/ M
| 
