 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） . s) k6 M' Q+ e- C$ M; X; d; N
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
; @. J& q$ f) q, F+ N! L+ D 
9 k, }+ w- w) j& B3 ?, w) B
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。  : ^- m* g0 G7 O3 Z3 W- G
安装nppr包
6 q1 M: @0 c( | v3 y 可在github（https://github.com/chaoxv/nppr）获取nppr包。
0 W3 m) v' f d, j5 E* [) _! ] #下载nppr包 , B: O* B9 x$ s- {% Y1 a
#install.packages(remotes)
" i) b4 N0 V4 h+ O' { remotes::install_github(chaoxv/nppr) 9 N: m$ O0 n$ b$ H- |' V6 E: C
#加载相关R包 , N6 C" x6 I( J6 [7 h5 }/ P
library(nppr) 7 \/ h+ ]5 \3 \2 M3 t! c
library(RCurl)
$ A4 D; W/ e1 M# h library(XML)
0 h; W5 T. r G" ]) J' w library(R.utils)
; {5 W1 f% B, b- I4 ?8 u library(tidyverse) 9 ]1 Z& a, Q6 y, U1 |
library(lubridate)
: Q+ p# R( t3 F( H! r% A# w- b 使用nppr包下载海洋遥感数据 6 W) s; U9 Q" T3 x* c+ K/ P) @3 w
nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。
1 f1 g( z p3 w 
' x7 ?, N. U0 @! h' \2 c 以下以获取海洋NPP遥感数据为例作个演示。
. e* e8 f& o* y% ^2 u #创建工作路径
* o' T- x( Z0 y U# i8 E yourfolder <- F:/R/nppr/vgpm
# p! d) N4 C" A* p dir.create(yourfolder)
9 D, V; Z: H* @6 x2 K #以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm
2 j- @# x( C* g get_npp_vgpm(
k; y6 J* v7 A" r s file.path = yourfolder,
* N v. x" T! d# a grid.size = low, #指定low或high可更改空间分辨率 5 M3 m0 @3 k8 a0 }
time.span = monthly, #monthly代表月平均，dayly代表8天平均
2 p4 }! ?7 s. w satellite = MODIS, #选择卫星 8 p. }7 C. U' z
mindate = 2016-01-15, #指定时间范围以下载遥感 * U; V. _4 P/ A, r. A6 T
maxdate = 2016-03-15
( N: K9 z+ V" q- g4 v )
9 c6 o* G$ G: ? 
/ j- a! Q- Y- j" w 在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
' J$ G {8 p# r+ I 使用nppr包进行遥感数据格式转换 9 _+ p" E1 l# S. Z" Q8 V( o4 E
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 ' [# y: @/ ?. z6 t( ?
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换 C9 v* p" e1 _: D( F5 h
yourfile <- paste0(yourfolder, /201603.hdf) 0 `( k& J8 _4 h8 i; I" j8 T
vgpm <- read_hdf(file.path = yourfile)
. S f6 @( _, Y) z0 R head(vgpm)
& C% k) G, L+ S) Q2 ]+ [& c1 d write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
1 D' ?+ b+ P8 C% c& K Q 
. W% I/ e2 c, B 转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
3 X5 G! A$ r$ s8 w 使用nppr包匹配目标经纬度的遥感数据 L, Y. k0 y1 m' b' r* \
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 , h1 @7 r; \) H; k
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP 7 i( G. a7 R& D6 g
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
- J9 @- Y% {% `2 e6 c( j+ z5 p #或者同时指定多个数据，不再多说
0 D' n3 i, j% _% E' Z; l! c mydat <- data.frame(
0 k x" }' n9 J- x- i lon = c(120, 112, 116),
j8 O! u2 a( {3 E' e lat = c(17, 15, 18), 8 g2 w P6 _+ f- V% i* f
date = c(2016-03-04, 2016-03-07, 2016-02-04)
; u& X% k1 T8 z# b1 o ) 9 Q1 v) z" W- P0 j, h5 j
match_df(mydat, file.path = yourfolder) 
i \& ]. h: }3 ]- C 绘制遥感地图 6 {. ]+ f% f6 N" n7 f/ ?
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。 : g3 x! x, t. g) b, S' M) }' E) p
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 ! w- v- G+ d- Z; G6 M6 c: o
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
5 ?7 K0 K N. F0 k( e library(viridis) ! P5 C$ y z; T0 g" \( G2 V
library(ggplot2)
^; Q" p) e7 e( p- X+ v' [. E' N ggplot()+
9 U8 ?+ @$ ]8 s4 p9 r- [ geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
0 M# @% `$ Z" ]$ o2 ^' ~% U scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ 6 b m% J& Q, p( c, l
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
% Y: o/ ^2 z: m0 ]0 ~6 X! p 
; }$ m5 @, l' n* d2 [. \1 M6 G ] 根据时间和经纬度列表匹配遥感数据的批处理
; R v: b5 w* T' S 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） ' Y5 x' m/ u# {1 k7 t% ~$ Y
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 $ H- A7 v+ a. t; u/ L! \
- k, I4 b2 D. J$ s
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
, r' X! _2 {* A4 @2 f- }: n ##如下以匹配SST数据为例做个演示
" w6 w3 I! I: M. `8 e dat <- read.delim(data.txt)
' b% i9 e: x- S# z4 K Date <- unique(dat$Date) #获取日期 8 r5 s+ W6 f4 U* W
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据
1 S) g% @+ n4 d( @8 z' U: h dir.create(yourfolder) - G. }9 C8 E; h' k8 Y7 S/ q! p
#通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例） ( b1 d' D9 a/ D% H: \& J: D+ j. v
for (i in Date) { 3 D' s# E0 [' D' f, o: x8 q' u
yourfolder <- paste(getwd(), SST, i, sep = /)
2 J9 ]5 h' @; B1 `2 X4 C dir.create(yourfolder)
0 u& Q4 R6 W% q) t get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) 7 S: e0 e+ m' U) w6 K
yourfile <- dir(yourfolder) " J- a: g4 x5 @* {) X/ E* F
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) & D$ K- G# P% e5 n$ W7 h
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
8 |, g% c8 C7 E& r }
8 x5 ~8 c3 B0 e #再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） % [/ J0 ^; C$ ]
for (i in 1:nrow(dat)) { , B0 R; y9 A8 I& ]1 m t. F
Date <- dat[i,Date]
% ~5 F* K* |; Z yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) 4 j9 H' v4 b W% i' {0 o' J
hdf <- read.delim(yourfile)
" b. n# {/ }7 Q6 \: |( N 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), ]
+ P2 _3 K( Y9 v5 k- h) _* i 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), ]
( ?& \& n. X' c5 J dat[i,SST] <- mean(hdf$var)
( Z4 Y Q2 H |+ A9 m } / K r& i$ J' u; D+ S
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) 6 d, \* `+ b J! g+ g
9 l4 a. l* L4 H& b" ^3 z
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。  ) W$ R; n5 z! S$ I" K1 x4 E
/ J0 s' C7 z+ D" g7 n
( K" q ~& S% X8 U8 V5 n5 E# G1 U# m; B- M+ ^
; q& X. C2 h, E; t$ R
| 
