 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP）
, v' D- d7 @: \" O' d$ y* M Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。 % v) |' v2 {% `" q! H
3 X( J* }" D+ F/ P4 w8 M5 t
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。  - x) s3 ~: U2 K5 H- Q b: K
安装nppr包
9 q u4 ]* X$ p7 n/ s' p 可在github（https://github.com/chaoxv/nppr）获取nppr包。
7 I& E' K- R: Z0 q1 @' B1 _ #下载nppr包
0 O" A% P* U. `2 M1 h #install.packages(remotes) ^+ b+ X. I& _9 U- g2 [$ ?) H
remotes::install_github(chaoxv/nppr)
3 I7 X. e& e. T' S$ [ [ #加载相关R包
' y3 F+ v" D, w9 e library(nppr)
" g0 M: x/ z7 p) w2 O6 O library(RCurl) ! ~& n( X! x" |+ u+ S$ e! i+ P: p
library(XML) 9 A+ r. ?9 |' N0 x9 r4 c6 F
library(R.utils)
: d% Y# G: c+ f; ]$ E library(tidyverse)
" ^5 m* J. @" P library(lubridate) 5 s1 S! O& A1 _1 G( J6 I- Z
使用nppr包下载海洋遥感数据
5 z. |+ P- Z! }0 d nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 2 g3 R2 g ]; v4 B' Z$ P* W) {
, I. B- |) t, T* w2 K8 x 以下以获取海洋NPP遥感数据为例作个演示。 & m h Z1 T, W/ J) @3 K
#创建工作路径
" a6 X' `3 W8 ]& g$ Z yourfolder <- F:/R/nppr/vgpm
6 t+ l) }4 [! F dir.create(yourfolder) & j) y6 P2 o$ X' Z1 E& {. `9 s7 b, o
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm " w% a" t% v/ ~- b6 N* F
get_npp_vgpm( % k6 y# s6 B" l5 J
file.path = yourfolder,
# r5 |8 _' C0 s- @3 u9 S% a grid.size = low, #指定low或high可更改空间分辨率
- a7 q. D) f, c, f1 C4 S time.span = monthly, #monthly代表月平均，dayly代表8天平均
5 M8 d- e5 o, R satellite = MODIS, #选择卫星 3 u$ M0 |& r0 z2 Q. Q+ H1 l1 {
mindate = 2016-01-15, #指定时间范围以下载遥感
) W$ S, P' @7 l( E+ n0 n9 i j maxdate = 2016-03-15
- d6 C* v! D) I- z )
- s1 T; ?0 ^( U) W5 u1 b 
) Z2 f* o3 q# t& p. v 在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 4 E5 m( W+ r- q6 e' i i2 I2 i4 |( y
使用nppr包进行遥感数据格式转换
$ Z# f; z9 H+ e# u+ f 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 " z" M6 o7 b0 P( h% e
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
. W# s. l, d& a/ W7 p. r yourfile <- paste0(yourfolder, /201603.hdf) & D: T; D" f/ z
vgpm <- read_hdf(file.path = yourfile)
# X1 Z5 A6 X9 i6 Q" m head(vgpm) : A$ d( Q$ W4 {; O
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) 8 ]/ W2 g! Q# m6 l
) Q6 J6 a0 ]$ z A8 f/ g3 t9 a 转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
X* J4 O8 s8 K) y; `( E2 J L 使用nppr包匹配目标经纬度的遥感数据 0 U% ]( x& R6 x7 ]& |
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。
Q# ~( u, g- k, k2 P. M# f #获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
1 _! @# F+ x+ X C! e6 U1 o3 t match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) # o5 H7 I9 e5 v2 \; P
#或者同时指定多个数据，不再多说 ' B5 K9 |0 R' h- n) J3 L
mydat <- data.frame(
% A; c) A- v8 `6 M lon = c(120, 112, 116),
# `: V3 G2 @) u8 C; N8 r7 ]4 P lat = c(17, 15, 18), ; \2 J2 D& C( r- f7 M5 v
date = c(2016-03-04, 2016-03-07, 2016-02-04) & h; ?* U6 K B( M8 ]
)
$ P8 r/ e( w3 ?' `+ k4 \& H" s; s match_df(mydat, file.path = yourfolder)  " T/ d$ f4 A' @4 z( \
绘制遥感地图 / @9 [, F: _; {) e I7 v d
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。 ! I& V, w1 ?& B; ^. a1 N
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 ( y/ Y& y q" q! W7 p! z) S
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
& R/ W5 a0 c, y3 L. K/ @0 S5 G library(viridis) ; e/ o% T# v& ~" d9 C
library(ggplot2) : G0 B9 k! ?3 C3 x7 z
ggplot()+ M& i7 Q+ |# v; P4 I
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
& t. w# o$ ^6 A4 T scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
$ \7 x1 Y- q& O, ] labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
0 W- B! H0 P7 _2 j 
& Y" @6 `5 _' Q P
根据时间和经纬度列表匹配遥感数据的批处理 6 r1 b* z4 z% s7 o; l6 q" H3 d
实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
9 t9 L r1 Y' l7 Z; I' O$ t 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。
! `% Q! L8 Y* w+ c! H! }1 H% ` 
" h5 Q' O0 p8 b. ^+ e5 } 随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。 " v3 M% F6 R4 [. F3 F. j) X9 y# T
##如下以匹配SST数据为例做个演示 - i A; ^+ c6 m6 x. P2 A8 l$ i8 ~
dat <- read.delim(data.txt)
" t5 Z R) q/ ^: N* { Date <- unique(dat$Date) #获取日期 $ J% w' b2 o: X6 K N q) R# p
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据
8 N( O) e2 K7 [3 L& \ dir.create(yourfolder)
1 V+ e/ @7 _: [ #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
, J! y2 ]0 v$ _/ H for (i in Date) {
$ R" U L, R6 e$ H6 P/ v6 w, U yourfolder <- paste(getwd(), SST, i, sep = /)
- r+ u# B# e3 N dir.create(yourfolder) 6 Q) \$ y8 b) F
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) 2 d; J5 m- Z/ o5 P0 o9 t. l& r
yourfile <- dir(yourfolder)
4 _7 C( n- D" S l$ f' M hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) % Z( b! ]0 |# m
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
) W/ \. i0 Y! r6 I6 A/ d6 f } . ~4 M4 C. N4 d5 }7 l5 a" F+ Z
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） + b! s: \, V. a3 q. D; s% I3 ?- p
for (i in 1:nrow(dat)) { $ s, {9 t" |9 q* D& G& I$ T0 V/ j
Date <- dat[i,Date]
& {* b3 f7 O3 e2 u: X# Q3 r yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )
. [# P( y2 I$ ^/ L# } hdf <- read.delim(yourfile)
% \& U) a( C- V9 ~! 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), ]
5 Z5 u/ v) t2 r- w7 p2 O8 r 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), ] : d% A% o3 i0 n" Y
dat[i,SST] <- mean(hdf$var) 0 V; a6 [, ~4 r
}
l% `$ z& y" q" M write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) * C: m, F$ @; x" J8 y2 `# | w
- J& F q# `3 G% r8 p: \7 r
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。  , Q4 l+ E( U4 Q/ t n, |& i
' b1 e/ @, s1 o; e8 g' q# R3 f$ F9 X
# @ {6 A8 @3 K7 E _( R4 R4 H: P8 L' J4 D/ @
& W: Z$ z8 U. w7 q
| 
