 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） * a" P- Y" K+ [* h7 Z- w
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。 4 b) L$ v, p, D3 G' J
: G) k5 E6 y; e" V
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
; f! r8 S n8 ] 安装nppr包
! f9 T# a5 X! ^8 h/ ]& @ 可在github（https://github.com/chaoxv/nppr）获取nppr包。
: W/ G9 | [# y$ K* `: k9 `# _ #下载nppr包 7 r, Z% N. v6 z6 `
#install.packages(remotes) & _. m7 K5 b5 G. F! i3 V' }
remotes::install_github(chaoxv/nppr)
' m/ v' j$ g* I #加载相关R包 ) W3 m' S% b% C i' f
library(nppr)
8 q Y+ [* h6 G, h9 I+ `# {, d library(RCurl)
9 a0 Z6 B- g! y4 T- _* x5 {- c library(XML) 2 y1 U& G& U6 ]9 I
library(R.utils) 4 g* R+ _& j M* ~ v- B
library(tidyverse) 7 J. b( w' _! P1 g$ r0 |0 ?
library(lubridate)
+ L. R8 e9 z5 H% @1 I' A 使用nppr包下载海洋遥感数据
/ c: ]% r' ?" D" k nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 , D: g5 I) }& A5 M8 r
/ F6 p( q: l# V- Y9 f0 a
以下以获取海洋NPP遥感数据为例作个演示。
1 y7 {$ e+ X1 a) t* @ #创建工作路径 3 D; v5 T' d! K3 n4 s
yourfolder <- F:/R/nppr/vgpm + @9 T2 v4 t2 }# M
dir.create(yourfolder) 7 z& c' J8 I \+ r* X; W" ~
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm 1 m0 O6 c. L! g4 Q
get_npp_vgpm(
/ ]4 S2 X6 R7 ~' A) v! v( l file.path = yourfolder,
) T! K* @( n9 J$ K w1 p, | grid.size = low, #指定low或high可更改空间分辨率 $ p, i3 Z1 }+ A
time.span = monthly, #monthly代表月平均，dayly代表8天平均
4 m/ c4 L/ @& [ satellite = MODIS, #选择卫星
% ~) z% o) M+ g: ^4 [ mindate = 2016-01-15, #指定时间范围以下载遥感
r: J9 H2 L0 M/ ^ f5 |+ b maxdate = 2016-03-15 ) _/ \# ]- T) R
)
, K( l# u# m; q 
! Y* T9 Y6 Q3 l1 W* k
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
% g8 b9 G* ~' c4 } \1 ^% @ 使用nppr包进行遥感数据格式转换
& }1 V- M8 Q/ { l- G E+ K 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。
# L3 z7 f8 H" M #将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
1 |8 g" ?8 Z5 h9 m* P' \ yourfile <- paste0(yourfolder, /201603.hdf)
* d% p$ Y+ Y& U' U. ^ vgpm <- read_hdf(file.path = yourfile) $ x% w6 Q$ X: H9 c, c( z" Y
head(vgpm) ; q1 m9 W+ w. J# U: E7 p' b' `
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
2 l" X0 A/ E/ A5 J4 N/ w7 `" l$ U 
0 m0 h8 [4 d- `! a
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
0 G2 E. R* K# M& ~! s 使用nppr包匹配目标经纬度的遥感数据
, L5 x5 z. J; H! U1 w 默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。
5 S, D" [8 C S( C, ? #获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP 0 W c' Z) c! v0 m5 r, f- J
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) " f' O) I, t8 g v9 D# s
#或者同时指定多个数据，不再多说
q: Z% ^# s/ k! r/ R& r8 N0 m# W mydat <- data.frame(
3 C0 k5 e+ F3 S6 g& k lon = c(120, 112, 116), , ]8 \: A) Q1 ?& ^+ K9 x/ s
lat = c(17, 15, 18),
. M' X" T( o- ]! b: F) a1 b date = c(2016-03-04, 2016-03-07, 2016-02-04)
' G" S+ F5 V$ y% C* e ) 8 V- r7 N8 v M O, U
match_df(mydat, file.path = yourfolder)  ' I" \6 i* `5 f# K, X* z
绘制遥感地图 0 F+ i$ N# |/ b
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。 + a" Q: w g3 o% W1 l) n
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
' y4 d% L v" M' w- O1 ` #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
" D8 r, K5 Z+ O1 H library(viridis)
2 m7 B9 b' X1 p library(ggplot2) , e K. V7 q; O% ?8 w1 \0 H
ggplot()+
9 Z" \; O8 T% P9 T( P geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ $ r, m0 ~; m7 J& `1 ^' N- O
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ ; Z% C% E5 j! H( W1 A
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) " J! u+ ]) z3 G8 C" x
/ T/ h" x6 F, R+ ~, n. x4 Y& w
根据时间和经纬度列表匹配遥感数据的批处理 6 c/ I/ A4 S' S" Z
实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
, N% v) ]/ J# K" z) E+ p$ |& `3 ]' ^ 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 2 i A9 Q9 S1 m4 o+ e/ x K9 x3 [
& q: v% w. j, a* `0 N
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
. K# F2 a. P% J* S7 Y ##如下以匹配SST数据为例做个演示
6 j1 W" ]9 t# e5 X: n* N dat <- read.delim(data.txt) 4 P( J' c; V3 r9 T+ L
Date <- unique(dat$Date) #获取日期
?0 z- b! |6 n6 `/ E yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 5 C% O8 P+ Z( F# ?% z8 _. E. C8 {6 H
dir.create(yourfolder)
. r: ~5 I3 ?; Y6 A M- X #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例） - |9 x& H4 c3 y, \+ y
for (i in Date) { * I' B# F- x/ _4 F% ^ O3 P" u
yourfolder <- paste(getwd(), SST, i, sep = /)
+ j' N3 U" I7 D& W dir.create(yourfolder) - Q7 P% q, W% u% c
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) - W5 T2 Z: {# C. j
yourfile <- dir(yourfolder) , m! ]: e1 [- S( ?6 v" E
hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) ' `; a Y& B6 Q. F
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) - c$ z8 P# T/ K3 M& s4 B% s
}
" ]- A2 D7 Z5 o9 ~8 d& _ #再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） # H2 X5 o8 K8 Z
for (i in 1:nrow(dat)) { " J5 S9 {9 u c
Date <- dat[i,Date]
: n6 M6 E( {* E! ] w ^ yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) * Q/ c$ P! F) n5 g* f7 D4 C8 q4 m/ ]
hdf <- read.delim(yourfile)
0 z8 M* A6 r( N4 e: o 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), ] ! g9 Q/ a# E' o
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. Q6 _& d9 G0 u7 |$ i
dat[i,SST] <- mean(hdf$var)
& A5 K: X! X2 e1 d& \5 s } " |& Q- ?& l* G) U& R
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
% y! z: l9 `" t& c& x, Z) V 
. M7 V, a4 O6 g) }9 f
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
1 Z, P1 t. f6 F/ c
$ j1 j8 W- Z6 K, h$ r% i% ]% O: z0 M3 \$ d/ D# A5 }3 R3 _
+ n5 \1 T, d8 i2 G
* G) ?/ C3 H5 [2 d | 
