 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） - ^& s$ \! k, G$ d7 X1 u
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。 f' M% _8 b6 c4 J' {. j% z3 p, i1 t5 @
, S" w$ U$ {3 D# u' G7 k f 本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
( o3 H9 e' Q6 K! }+ R7 x R8 k 安装nppr包 6 N. L% H6 W* k
可在github（https://github.com/chaoxv/nppr）获取nppr包。
5 _7 g9 a4 \# `; E/ `7 ~ #下载nppr包
( {2 X) C3 y9 K5 [ #install.packages(remotes) ' E O, Q8 W7 ]" j' R
remotes::install_github(chaoxv/nppr) ) ?) a* \3 m7 N# I7 h2 W
#加载相关R包 0 O& B' p1 N4 c2 U
library(nppr)
4 e! K! ~; x" J' T9 e/ S library(RCurl)
9 q5 T# H! X- A& G* i library(XML)
3 c. S& W6 x0 X0 ` library(R.utils) 9 }9 `8 {8 W2 o; Z7 J
library(tidyverse)
: h/ V: D! C8 }( S; q library(lubridate) 2 d4 \5 c, r! z
使用nppr包下载海洋遥感数据 ( a1 K5 U# ~# G. A# }! [
nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 / Y7 ?5 K x2 m
0 J- s& b) m% e4 h3 h% M 以下以获取海洋NPP遥感数据为例作个演示。
F+ n- p( X. l& F$ T( d+ O6 I #创建工作路径
6 ]+ q1 G! z1 d8 \1 [6 Y x9 | yourfolder <- F:/R/nppr/vgpm $ F% G. u6 |. F9 Q
dir.create(yourfolder) ( o2 U6 Y3 @& h. u+ Y+ x- l
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm
1 O# G6 B: c/ H I f5 S$ r G7 | get_npp_vgpm( & ^7 G6 e/ w# E* }+ H" s
file.path = yourfolder, . o% l" A2 Q% q
grid.size = low, #指定low或high可更改空间分辨率
2 Z1 ] z$ A; y6 g' M* X time.span = monthly, #monthly代表月平均，dayly代表8天平均
! F X4 K; T0 {! M9 ~5 f% d5 a satellite = MODIS, #选择卫星 ) U( c. m. ^" s0 ?4 ?5 ?( }
mindate = 2016-01-15, #指定时间范围以下载遥感 @8 _0 P! N8 w8 ]5 N9 s/ F
maxdate = 2016-03-15
6 y: ^7 X# E! A% d. ~0 _ )
+ D7 L# t; j! |, r6 ] 
3 L! d3 p2 M$ E @# ] h( { 在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
- |, Y' D: C- s7 l1 ` 使用nppr包进行遥感数据格式转换 3 n4 k2 @% {1 v- q: w7 J
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。
9 K3 `9 a7 j' `0 y6 {1 p c( {4 T #将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
8 A: x0 X0 A! }; D8 o yourfile <- paste0(yourfolder, /201603.hdf) 6 O. ^' u2 Z+ t9 S( g# H
vgpm <- read_hdf(file.path = yourfile) & M5 P2 u B! b. j
head(vgpm) ! O. q8 f0 ?6 L' B* e, v
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) * @! G& u | w; P- f0 K5 c9 _+ @
# N ~4 i( n8 u6 I& e! W% [
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。 : b1 {( s9 x- k e% o3 ]: l6 l
使用nppr包匹配目标经纬度的遥感数据
- J8 d+ Q6 g8 |# v9 T3 p 默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 ! g2 z' B3 ` j* p0 V; w. f
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
$ u: r; S; s [( l match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) , p# x( z9 r/ Q4 j9 B/ _
#或者同时指定多个数据，不再多说
5 F* r% i* q, w mydat <- data.frame( 8 J% W9 r' z' }
lon = c(120, 112, 116),
0 h4 I# o+ A$ q& S4 A8 E lat = c(17, 15, 18), " W# w* g! [. N" u; [7 u1 q
date = c(2016-03-04, 2016-03-07, 2016-02-04)
6 O: n; `) P0 [8 I ) 7 ^" Q4 f7 G: ] ~1 J& l
match_df(mydat, file.path = yourfolder) 
+ g; w" A; k) j1 \' K# z+ m 绘制遥感地图 " x, q6 B$ K& R: T' e
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
" s1 H5 ?1 u; t9 L #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
$ t( R* |! x& ?8 R2 \- b. c #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
' a, Z, r+ H, E/ v# F; | library(viridis) ' R1 B5 w. N+ ]" N) d+ B
library(ggplot2)
/ K) l: \+ S6 A; b9 x4 w ggplot()+
( Q- b8 T4 ~0 [& h geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
9 N8 s' |& J) {; N# ~/ @ scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
& v6 p; C$ _8 R% E' @ @ labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
! j8 g' W* ^& S- |9 N 
" c1 }4 f8 b( q+ E% A L- E7 @ 根据时间和经纬度列表匹配遥感数据的批处理 ; u2 y' T0 H: p2 K, j! N+ t, o
实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
( ~: x2 g% f; ~; M: b. L. Y0 O( T 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 / m9 b5 q* D; C$ ?
. E; n; m3 F9 Q' U5 J/ g/ n
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。 + U% _6 [* }# D6 l
##如下以匹配SST数据为例做个演示
' i0 l7 _2 d8 H7 M5 a& ]5 c dat <- read.delim(data.txt) % v) t& O: \/ F6 ]+ [* C( N; r' y
Date <- unique(dat$Date) #获取日期
! y6 H/ [# k3 `' M7 N* w yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 , h4 e/ v& I7 Q0 S3 G! B8 p/ H
dir.create(yourfolder) & Y0 q$ r+ Y8 d1 O1 K( n
#通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
- C- \% L" n5 P D for (i in Date) {
" `" m# e: ?; f yourfolder <- paste(getwd(), SST, i, sep = /) * q( ], }* P+ \: M9 ?
dir.create(yourfolder) . O4 g4 X) k( V$ w
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
% |& j9 T, k% X6 ?5 `7 | yourfile <- dir(yourfolder) 6 f$ x+ r) k1 h" p* d9 o7 q
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
6 n' R( v: L$ Z: I5 G, k' z write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
( ?" h& T( I+ x8 @ } - H) x( X* z: f8 ^
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） : \1 J" N7 X+ v0 i$ \1 D
for (i in 1:nrow(dat)) {
1 h! K+ H$ T+ T Date <- dat[i,Date]
# x/ m0 c+ T, @9 f# y6 }; g- j yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) $ y9 M( Y: j3 q7 A
hdf <- read.delim(yourfile) ' r& B4 N, ?4 X t T' b
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), ] ; D2 R: e* D3 t& g
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), ]
) O0 k& K. {$ d4 d dat[i,SST] <- mean(hdf$var) 8 f; u V+ N" B
} & A( e* T2 G2 g3 I2 r; A
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) 4 N1 F! t+ @& w& W- g9 V
" a% {5 V+ C' m- F% l6 ^1 b' I 输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
* b4 N: D: n+ u. D4 ]) k$ V" E1 W' L& H9 h# t, G
! ^# n1 x9 T+ m- r( `6 C
`. V$ k! o8 R; E$ p5 M/ o4 n: P7 v3 }( Y! C
| 
