 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP）
& X' Q5 f; }. [" R" }( X Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
# e& ` G9 V: @ 
+ A# B1 F6 E, a
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
E& W# u1 A% f- w8 ^1 }" [ 安装nppr包 : s0 F' v6 N; c9 F5 V- g) Y
可在github（https://github.com/chaoxv/nppr）获取nppr包。
4 n& b- k7 @% C8 @2 I0 Y4 K #下载nppr包 ; H# C# a# G- k
#install.packages(remotes)
, h u. y6 D e* T remotes::install_github(chaoxv/nppr)
. R; E" ^. H- h6 _ #加载相关R包
0 I. M2 A/ b* Y; n' p library(nppr) ( D$ k$ I, p7 s: _3 h
library(RCurl) 7 S% q6 h5 e& J
library(XML)
4 ]& Y, w4 v. K6 R# L. x library(R.utils)
7 E, k. Y4 o# S: A8 o( s. D library(tidyverse)
9 g( S! \0 e' I7 P3 V% | library(lubridate)
5 E- T; ?/ a( l# @+ k' U% E, m 使用nppr包下载海洋遥感数据
5 n1 x( p) ~% h0 ^6 m/ s nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 k6 V6 H6 N: e% j' r
) D" a5 Z) C0 ?, ^
以下以获取海洋NPP遥感数据为例作个演示。
( I+ i+ Z0 ?) ~1 H4 I #创建工作路径 1 L1 f" q* z& l! c. _
yourfolder <- F:/R/nppr/vgpm
3 N$ `( y# N9 c, B& r' |3 @ dir.create(yourfolder)
& k$ C, u. d$ Y; F* M5 F #以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm 2 \; G& \% ?- c7 l
get_npp_vgpm( + E I# t; n0 W3 Y
file.path = yourfolder,
9 J2 [* \2 A# Q7 t$ w grid.size = low, #指定low或high可更改空间分辨率 . v2 o/ Z) R q! I* U5 G
time.span = monthly, #monthly代表月平均，dayly代表8天平均
3 s7 u/ H+ U; w0 U5 O6 z/ F1 E7 f# M satellite = MODIS, #选择卫星
; T; r7 K& Y+ n0 N/ Z mindate = 2016-01-15, #指定时间范围以下载遥感 - V- w7 c0 ?" r; t. t5 F3 b
maxdate = 2016-03-15 ( P. i" y0 P; ]4 x
)
' _4 r/ t% C7 R. k$ x 
/ I" t2 M2 X9 h/ S T 在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
( ~6 i3 m5 {8 M1 I 使用nppr包进行遥感数据格式转换 ) Z8 Z7 B7 H* u4 l. q1 ?( ?" r
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 ' Z" D. D/ Z+ m
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换 0 ]7 ?& }4 d& B5 ]
yourfile <- paste0(yourfolder, /201603.hdf)
' d% A, [6 [: }- G- ` vgpm <- read_hdf(file.path = yourfile)
: M! P9 Q& ?# e' W0 d- W. c+ s head(vgpm)
" r M% w: b, a( e& L1 C% ?6 ] write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
6 u1 ?$ N: F0 S+ C2 K& [1 [ 
# |! t+ s: Z" p5 q1 F0 P
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。 + x: |% p# v d( @
使用nppr包匹配目标经纬度的遥感数据 " E0 q5 h+ _0 U' @
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。
# L1 _" c+ \1 \$ n; z. Q #获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP 8 r' T! @# c3 \7 j/ \
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) / Q7 M; z6 Q$ Z9 s/ W
#或者同时指定多个数据，不再多说
?" r, K& l/ p& b# G8 ^ S' k mydat <- data.frame(
# M8 b( c5 n1 M lon = c(120, 112, 116),
$ h; N" V- q- _ lat = c(17, 15, 18),
3 h6 O8 v- x( I; A# a/ w# i9 t/ V date = c(2016-03-04, 2016-03-07, 2016-02-04)
" z2 ?3 z" B0 m% K I ) 6 Q2 m$ B4 U$ j* i$ ?" {$ ^
match_df(mydat, file.path = yourfolder)  : \) P [0 k8 g V5 T3 V( m5 Q0 K6 Y
绘制遥感地图
4 Y: u! }* e- m+ I: V4 Q' E6 J nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
$ K7 @7 Y1 C% J L& h #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
# W; _8 ]: @$ I2 @ #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP . \; i. a1 }' X5 P0 {
library(viridis)
n$ u f9 b1 K* a7 m library(ggplot2) % {4 ]; p2 s2 b2 Q$ x+ N
ggplot()+ : ]# A( {6 e3 i& w$ h. M
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ 6 m3 u/ G. q" O
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
) j' M, Q' \' J labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) ! w# v4 S$ `1 C0 j0 S% t1 ~2 z
+ S7 o+ V7 S' s
根据时间和经纬度列表匹配遥感数据的批处理
6 m; M' h, B# N" ` 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） ) k( L2 a D1 C: G! ~% X" y8 @
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。
( ?: \- j) P- N* o4 F 
; ~9 \2 R5 l/ F; p V' j9 I2 X
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
6 L. U/ \; z2 f" g5 w ##如下以匹配SST数据为例做个演示
1 K4 D! E3 y1 A _1 ]2 f1 u dat <- read.delim(data.txt)
/ s6 e* N% o2 a$ _+ z. {. M7 U, [ Date <- unique(dat$Date) #获取日期
# |$ F0 E4 J' A* {" {6 L yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 ' D- }$ Y5 l+ R
dir.create(yourfolder)
7 ]5 [% M: r/ ?: v #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例） % [; u9 e, _; L) T8 k2 `( q+ D
for (i in Date) {
: B+ ^% F, g; b! Z2 _6 ]$ T yourfolder <- paste(getwd(), SST, i, sep = /)
3 L# p3 Z7 U2 e4 O6 s9 d% J dir.create(yourfolder) 1 R+ `# x. h& r3 J* ^
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) * R) r/ l2 n1 `+ L T
yourfile <- dir(yourfolder) $ r% Y4 c8 A! v+ c& o) u
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
5 I: |& p+ O) B9 O1 ?5 R write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
9 d' `4 c3 z# A } / S( Q' A4 F* I5 k% \4 l
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） ' e& F( B7 u! `- e/ w. z
for (i in 1:nrow(dat)) { 9 S2 A% @/ M/ y) h0 _
Date <- dat[i,Date]
. ]3 c2 f5 b0 a* s% d yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) ) }" K( o4 @" H/ ?6 e
hdf <- read.delim(yourfile)
. `5 I$ ] V2 H$ D- Z 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), ] " X% V6 m. b! V' E) S% s
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), ] / r R% [- i& I" d2 U
dat[i,SST] <- mean(hdf$var)
& V' ?7 q- ]; g } , s* f: l' N# W% a9 @4 R
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) 2 r q. I- z$ L$ o! Q- Q' ~ ^
p" F, u( d; q0 y 输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
6 b0 w q4 A7 u% z
; H: `& H$ M- M4 j4 v# Q: e8 U: H! h2 q' Z4 E2 v; t
5 `/ }+ _9 x) A8 _) W0 Z
! ]: r ?3 H+ P: h
| 
