 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） G$ {; F' Y+ X, P; T
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。 : P u7 \; r! f1 r+ o Z
8 k' U5 G" Z2 `5 ?" z 本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
Y9 b' p8 y N* F# b 安装nppr包
. C6 f" B e r* I$ ~) E/ @ 可在github（https://github.com/chaoxv/nppr）获取nppr包。
* Y+ J) ?3 z7 u+ {0 J #下载nppr包 ; b# i2 r% V, K$ G
#install.packages(remotes) 3 c/ q! E# O9 i3 y0 m' m, K/ y9 N
remotes::install_github(chaoxv/nppr)
! l3 J% _2 B7 ? P6 k #加载相关R包 & |1 A- i0 O! j; q7 k
library(nppr)
- S. C* k, H& y2 T5 _& W. g library(RCurl)
3 ]# R$ Z, _! U% ?' I4 ~: n library(XML)
5 A6 B3 }% X( G; h library(R.utils)
+ N5 n- b$ G, ~' S9 X2 ?# r library(tidyverse) ( ^4 p9 v( U6 l2 E P, n& V
library(lubridate) " P0 W E& W( l% Y, i6 g( s) i j
使用nppr包下载海洋遥感数据
# _$ E* }/ h; y* k7 x nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。
/ J6 t: b9 A* U, { 
2 U( H, q! ~! X
以下以获取海洋NPP遥感数据为例作个演示。 4 e0 x; j7 ^3 _
#创建工作路径 L2 j9 G8 P- l+ ~& \2 \
yourfolder <- F:/R/nppr/vgpm 0 s" `* ?, l c2 V. A6 {9 I; n d
dir.create(yourfolder) - d: B4 V/ t5 r! n7 R- \/ Q
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm ! f3 \( v% M! }. K( p. {5 |: W
get_npp_vgpm( 1 {4 G+ z' W, U6 r! o' e
file.path = yourfolder, 1 w: P+ J |4 O: D: L6 v2 K% }* U
grid.size = low, #指定low或high可更改空间分辨率 3 }9 L2 a+ @6 e4 P
time.span = monthly, #monthly代表月平均，dayly代表8天平均 0 o' A/ Y4 t7 c# L* ~; }$ P
satellite = MODIS, #选择卫星 - a2 T( n; I+ g, X3 G8 ?
mindate = 2016-01-15, #指定时间范围以下载遥感
, l8 a2 ?! d" M+ Z& c maxdate = 2016-03-15
4 g0 ?* M4 H) e0 w/ {6 ^& k/ q ) 7 r' N7 A$ P- ?) T; B Q6 V9 Y
2 b+ {3 \% G6 D+ {$ e
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。
- E5 I3 I+ L, [& p" E b0 u7 T 使用nppr包进行遥感数据格式转换 # Q2 i" Q8 T% O- H: o4 i
如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。
: a7 F! P5 R: X; \# h2 G6 K- @ #将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
# T5 A# y1 u9 a8 ?) P% a1 ?# P yourfile <- paste0(yourfolder, /201603.hdf) 9 ~# O. e. V5 Z! B/ A
vgpm <- read_hdf(file.path = yourfile) ! v; u4 E4 l' d: S& c8 x/ v
head(vgpm) ) h8 _( `' S2 x/ ~3 y
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) # t: E' j3 k$ \1 ]3 O* i& c* X
# L- z& `4 \) @" [$ M# c5 t8 Q' ]
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
+ p7 v6 u- b" D1 k4 o 使用nppr包匹配目标经纬度的遥感数据 ) w* Q( t, ?* j8 Y0 Y* a" s/ w
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 5 I+ u' ^2 H! c: O+ G* @' {' B: }
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
4 R- l3 A8 @) X, f; b$ d0 y- g match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) 2 v" T. o! ?" E% w @
#或者同时指定多个数据，不再多说 + R2 I+ Q- [$ g% Q
mydat <- data.frame(
; P j6 _4 W4 y+ f7 e y! p lon = c(120, 112, 116), $ |( M8 l4 U) O
lat = c(17, 15, 18),
' h) R+ u6 `# [% \4 S- w# K date = c(2016-03-04, 2016-03-07, 2016-02-04)
+ A" |+ F; c* Z/ U+ w, E )
# x" ^3 K* W6 Q1 ~1 U match_df(mydat, file.path = yourfolder)  4 U$ @9 R5 v0 W+ z% ~
绘制遥感地图 ! q4 u+ ~ [! r
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
% [7 k# \5 `6 r8 l/ J #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 ) }# D( Y' o' a# e
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
+ V3 K& G* N( {6 S( p library(viridis)
% I! Q0 l; I$ E library(ggplot2)
. K$ [8 Q& g. u+ n0 ` ggplot()+ 3 |: `& o2 B+ K5 b% L. Q9 b1 ^
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ / j% z& [/ v( X$ _- y: j9 B- E
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
( _2 W) K! j2 S9 q6 x U8 N. r% ^ labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) . o5 _3 W4 U% K) J% H1 W
: T4 s; W% L; |5 e1 T6 ]2 X0 K 根据时间和经纬度列表匹配遥感数据的批处理
" q2 A1 @: w0 A2 |# M 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） - Y) g+ ^/ \. Z) @
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 # H* o& w0 G# T2 v
n- a* i' ^2 H9 P% c. |
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
9 t. U7 b) [4 Z3 Z ##如下以匹配SST数据为例做个演示 , D' x" [# j7 X& d2 m! B1 G
dat <- read.delim(data.txt) , C' l4 N |7 R8 v! W' M7 Q
Date <- unique(dat$Date) #获取日期 ' i4 e" }7 |6 e9 r, I1 E
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 % l2 q) z; G" _: Z! `: R9 I! I3 Y
dir.create(yourfolder)
$ U: E4 o; U& Q) @: L #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
7 b/ _0 Y) Y8 ?5 K2 d8 } for (i in Date) { & ~# {. _ n1 ]- w
yourfolder <- paste(getwd(), SST, i, sep = /) + N9 P7 l o R4 K8 {
dir.create(yourfolder) 2 _' F& o; ^' Q
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
3 `9 k' O/ [9 g, X( _5 E yourfile <- dir(yourfolder)
8 }7 y7 L' ?- a, M5 r hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
/ t! B9 V) B w& Q write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
: f( ]% s1 d# e } 5 @( f0 U) z9 d, ~
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） ' R( a4 _+ o) v3 j, e8 l! S8 k$ G
for (i in 1:nrow(dat)) { " m/ i% v2 s# H
Date <- dat[i,Date] 3 L- n" {) j, j9 [ I! j
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )
( X. w Z K( Y hdf <- read.delim(yourfile) - v' v5 e# ^% q0 a0 l: r1 L9 q
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), ] & R0 a# E+ q: R/ {) v. X
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), ]
3 `* n* b+ @' C) [1 l dat[i,SST] <- mean(hdf$var)
( Q2 Z! ~/ V( F- @! R5 F+ h1 e } 8 x0 X0 p/ V% Y% Q! Y
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
: A. |0 C6 U9 z' Q 
2 ]' i0 I* S$ P8 F) V4 h 输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。  8 `6 u' {8 M: S' K" a9 I. X* E; N
$ g0 I" d; t' I' ~" @# A2 }& f% g3 c3 B
7 y7 X& @! o, D- ~% ~4 @1 w9 I# j: z$ s
| 
