 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） / e( J' m: Z) |, m" e
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
8 p$ }) V! d7 D1 l 
5 J, }* L: }7 Z2 }- c9 ~ Z
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
0 L. S r& p$ ] 安装nppr包 * p) V' m0 B6 b" k- r* [/ D
可在github（https://github.com/chaoxv/nppr）获取nppr包。 0 M$ G3 R3 q2 Y3 k! V: F
#下载nppr包
+ I5 @* F6 d6 i0 Z& W; p #install.packages(remotes) , c! z$ V1 `# z2 N* G: x7 T) g
remotes::install_github(chaoxv/nppr)
: ?; s- M# j5 c$ l, b' T( r) _ #加载相关R包 ( a; R1 c$ B# k$ m
library(nppr)
1 q! J$ t; O1 J library(RCurl) 7 f2 b* X4 o2 h
library(XML) 0 v5 t% R3 d6 M7 b$ x' @
library(R.utils)
, s; _! |! S7 B9 _ library(tidyverse) : w, O& ?9 u4 {: R# W; b
library(lubridate)
6 R9 c1 j5 [, T' [, x7 r- l* L 使用nppr包下载海洋遥感数据 " i9 Y j/ n% R0 v4 F
nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 ' T3 F S m; @' C2 W
6 b m. J5 ]) _% T* ~& y 以下以获取海洋NPP遥感数据为例作个演示。
+ p1 Y$ h( D8 f) i. L5 R5 w #创建工作路径 - [* t+ s2 m6 X' A2 ^* @% A
yourfolder <- F:/R/nppr/vgpm
3 j4 ~+ r/ a/ N9 w dir.create(yourfolder) 2 |' O7 K+ M; w
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm 2 \; ?! l% ~' p% U: M0 W9 F
get_npp_vgpm( - W6 v% u8 O# q0 K; j7 I) {
file.path = yourfolder, 0 C! X& ]; w$ {' l7 V/ K
grid.size = low, #指定low或high可更改空间分辨率 " t @% |6 g: }. |
time.span = monthly, #monthly代表月平均，dayly代表8天平均 # y& d$ V+ X7 ^" U% t- U2 L
satellite = MODIS, #选择卫星
8 i2 z2 m) d- v6 T. y mindate = 2016-01-15, #指定时间范围以下载遥感
4 _0 m% A) e1 B9 I6 [ maxdate = 2016-03-15
) l$ p }, }- r ) $ _! v2 \* L! z/ Z) I4 i! X
% h8 o' M" P) I5 n6 i8 M! u( v: o8 M
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 & o9 `/ k. ?6 j& y3 Y9 B" S5 Z- {
使用nppr包进行遥感数据格式转换 0 s7 `! B+ P/ Y1 C* g3 g
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。
9 J0 R5 ?0 {- H; H+ W #将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换 F: g a( H$ k: v1 K
yourfile <- paste0(yourfolder, /201603.hdf)
+ U9 m0 \( F* i8 u# |# n/ k! C, ~ vgpm <- read_hdf(file.path = yourfile)
8 K. o k( ]5 E4 D head(vgpm) 3 I$ y. @" L4 `) g" Q6 `& l
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) $ `4 t6 ~- n6 p% g
0 m* B( }/ M5 [! L 转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。 1 p+ v8 W. j4 v1 I8 h
使用nppr包匹配目标经纬度的遥感数据 * b4 Q8 H4 k* `2 B- X' p
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 ) P# [1 P4 A6 F6 H- F0 _1 T
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP ) G# ?$ V4 E( G* N0 R
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) # n% |. S6 {9 ^9 D7 B, r$ v* x
#或者同时指定多个数据，不再多说 ) [ C: A, ?; n- {. h4 A5 ]
mydat <- data.frame( * |9 w+ \* v5 H, n+ n- ~
lon = c(120, 112, 116), ) H5 m/ J) [2 z% F4 m3 q; I. t
lat = c(17, 15, 18), 7 G: G% ?- |" n: v
date = c(2016-03-04, 2016-03-07, 2016-02-04)
5 d* S% q0 f' [/ J, D. k$ i )
% y+ E* F2 g& k5 r L match_df(mydat, file.path = yourfolder) 
, s5 s( m& m* J 绘制遥感地图 ' m# o; w- V) \- m4 v/ ]& ?
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。 1 H I9 f0 `* p6 t9 U
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
3 j5 a* @5 y7 n" Q& C' P# n+ A #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP 5 X. ~( Q5 E1 z9 R+ T/ g1 p
library(viridis)
8 R3 ?9 g% l& `# p library(ggplot2)
/ w$ ^2 C% |$ p& X) } ggplot()+ ; S! z( P: g% K1 |6 f2 ~
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
; n9 n7 {. P2 p. H- T. i" L scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ ! T0 C! ^0 f# i8 F
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
' U* b- ~3 d# e: O 
% ^# H, g/ s; w. g8 V& `) ^6 s; V 根据时间和经纬度列表匹配遥感数据的批处理
0 A" i2 Z9 z( N4 J 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
- v6 x3 z+ e7 E { 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。
2 E. ]' @& x5 \ 
. t( d, U+ J1 x& H; q 随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。 & m7 d1 t% n$ N+ n
##如下以匹配SST数据为例做个演示
( m5 C* _* Z$ X' X" X dat <- read.delim(data.txt) 5 U$ p4 ]6 ^% [. n
Date <- unique(dat$Date) #获取日期
9 I; z, `: L8 u. B yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据
' n. U- w$ N2 j3 z0 ` dir.create(yourfolder)
# m2 @9 E) a+ w% i" { #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
D: x; A2 e( R* ? for (i in Date) { : {2 e# J! N* q9 ]2 _ `
yourfolder <- paste(getwd(), SST, i, sep = /)
% i+ z3 e8 \+ ?* C3 i# G dir.create(yourfolder)
& }* j0 @1 V$ t/ d get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) 6 A5 h- b! c$ V3 l/ R7 d
yourfile <- dir(yourfolder) & T, @ P: Z9 R/ o
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) ) x! ?& K6 A- B/ P
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) + w2 |" H5 a8 n3 ?, v7 r6 N. h2 B$ Z
} - _4 @$ n, w1 R6 a0 Y
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） 8 F' a: t7 m8 R3 I2 d
for (i in 1:nrow(dat)) { ( r0 E3 x2 K& ^3 ~4 \
Date <- dat[i,Date]
" n) I9 T" m7 e# f# v9 V yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )
/ I0 Q, @; U @3 i; Z7 f hdf <- read.delim(yourfile)
6 w' P- Q- ]4 q1 s$ ^4 b9 D 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), ] 1 L1 @' ~; Q% Y4 e6 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), ]
K" h$ ]( B5 P8 B) P, x: \ dat[i,SST] <- mean(hdf$var)
0 l9 @- e3 d8 R& L% | } 8 l/ n3 A2 q% \* [
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
) u6 k6 b" E/ G0 C( D N7 H b 
2 a( X1 T& ^, w7 [* R6 }1 {8 ?
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
% U# T# y& v0 Y9 _; N
# _* Z. z: ~+ {+ ?1 N" u% y1 m# r: @8 b# X! p- @: m
( p% }- G r$ ]- _* A( S" D* A: M. _! L' K1 P; R7 l
| 
