 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP）
# U& g+ s. @$ i2 w Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
- U' Y$ d+ i4 }* \5 M 
$ P, y* D3 k; K: q 本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
) I# w. v7 E) a 安装nppr包
+ b6 p/ t# r m6 C2 k 可在github（https://github.com/chaoxv/nppr）获取nppr包。
6 L6 @8 _( ~, W( l, O% w7 C #下载nppr包
. x, Z! J# a: s1 v #install.packages(remotes)
{: j p8 P5 r' L, P7 ? remotes::install_github(chaoxv/nppr) 6 r0 v* x: c" Y
#加载相关R包
8 g9 A5 {( H. v& `+ P6 v library(nppr) * K: q5 K, ]( I# Q- Z! q
library(RCurl) # a+ c8 q4 f3 J5 Y$ l+ D
library(XML)
7 x9 I0 Y4 @9 d library(R.utils) " R5 Y: P* g2 y9 `) ]% n
library(tidyverse) - r2 X- X9 S8 O/ h# H
library(lubridate)
7 K c# S# o: K; p$ O8 v. j# N9 H 使用nppr包下载海洋遥感数据
1 e( ^/ u, G7 B4 Z1 @ nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。
, F: h! K& b& y# [6 W& l) o: M 
% Z/ T V7 ]! U 以下以获取海洋NPP遥感数据为例作个演示。 " t' A! u& ~8 N& B
#创建工作路径
. u. o! F# h' i0 ]1 ^ yourfolder <- F:/R/nppr/vgpm ) i! u5 p* N' B0 U* G
dir.create(yourfolder) / o2 Y$ {8 R/ Y' n2 C& r$ O- `
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm - { r' a2 { @: Z: _
get_npp_vgpm(
" W4 d* U1 y8 ?. e: F file.path = yourfolder,
, G3 ~5 O* c k grid.size = low, #指定low或high可更改空间分辨率 - q: _4 y9 C" s% d
time.span = monthly, #monthly代表月平均，dayly代表8天平均
- ?$ r: f7 ~* U3 a3 c; { satellite = MODIS, #选择卫星
/ s2 M. R' f' w7 @! K# z6 C8 V mindate = 2016-01-15, #指定时间范围以下载遥感
8 T# @* X7 e I8 ^' R0 w! a maxdate = 2016-03-15
& |4 B, X# L' n* {0 C ) 1 ]( R P2 ]% ~; l: W" w
$ N5 w& l- q7 e3 R2 ?3 [ 在这个示例中，我们使用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 b% }0 c9 a' I6 e; q) N4 ?
使用nppr包进行遥感数据格式转换
B6 Q5 o3 u8 d! r 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 4 `1 t) s" C" E9 X6 @* n% T( z
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
5 }# t- A' }5 c: Z0 I2 g yourfile <- paste0(yourfolder, /201603.hdf) 5 Y1 m* \3 i1 C
vgpm <- read_hdf(file.path = yourfile)
7 s* M+ d" j& x. X# z head(vgpm) 9 v) f5 G8 j" K5 d& u
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) 7 n' v5 c0 B6 T1 u. h1 H
5 T+ {8 L9 l7 x( J
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
3 h( ]8 Q; Y" H 使用nppr包匹配目标经纬度的遥感数据 ' K* p) q! E0 r# N
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 0 p, E# r# q& X" j& ^ z
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
) k2 H* l; w. L% r( I match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
q- W- s2 v- P #或者同时指定多个数据，不再多说
1 D: O' t5 z: f x& X% D mydat <- data.frame(
3 P1 Q9 B6 D/ z) a0 W9 ? lon = c(120, 112, 116), $ X* u" H8 I) ?* ` T
lat = c(17, 15, 18), 2 n! s1 L( r5 J- y5 n v
date = c(2016-03-04, 2016-03-07, 2016-02-04)
# G ~8 z! e) W/ `$ C ) 0 Z; F9 h5 E: w/ J/ p0 r
match_df(mydat, file.path = yourfolder) 
; O" N& m' S# O' A4 Y& d. {6 r0 c 绘制遥感地图
* P" a; }, |0 I6 \3 A8 }) q) }( ] nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
4 d& i* V& T* A# P& Y #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 ) c! q& a$ _- g: i( H% L# {% v
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
7 }/ R! c( D, G, x library(viridis)
1 v7 y1 Y) F/ o/ C0 k. S) `- G library(ggplot2)
& N8 U n) k7 u- I- z ggplot()+ : n+ u$ Q. O4 F2 d# |. G
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
: b1 f- j9 o) K scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ 4 Z+ o+ ~2 }$ U4 A
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) 4 s8 h& P `1 }: k- U
0 Q; e& N, q$ [2 K, `8 i% s
根据时间和经纬度列表匹配遥感数据的批处理
8 f# g! X" [1 J4 u7 V" B 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） 9 \& \/ t1 L2 J! C- s* Q+ \
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。
3 N* ]/ T, d+ A5 W 
+ x2 l* j K$ t2 e 随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
+ I$ _1 t" Q" C. @/ ^+ d ##如下以匹配SST数据为例做个演示 $ j' l* A* d* n: H" j/ n: V
dat <- read.delim(data.txt)
& ]1 W/ z+ Y; s+ U Date <- unique(dat$Date) #获取日期
: h# k. N" N4 r, f& d/ e7 k4 K yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 " K; ]2 p( F; V! R, J' J
dir.create(yourfolder)
4 ?. ^7 `' a6 O# p3 B0 |7 [ #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例） 7 |5 W8 n* F# g
for (i in Date) { 3 L) M `+ y1 {( [3 ^7 l
yourfolder <- paste(getwd(), SST, i, sep = /)
$ l! d% ?' Q% Y' u% P% X( S dir.create(yourfolder)
& s# [5 N6 d5 v$ ?1 }- x4 T get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
2 ?5 X: |* F% Z1 R yourfile <- dir(yourfolder)
1 O9 s' I+ E) r. q hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) 4 F, b) ^! @9 V
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) 6 R; N j) t H3 T$ i% X
} 3 P# e0 s/ w9 }2 R9 U: V* S
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均）
2 x% F+ h8 I/ V: e for (i in 1:nrow(dat)) {
: t4 m9 L) b {# b Date <- dat[i,Date] 4 \* a. o9 n2 K R
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) / R- v6 y4 f1 P8 D
hdf <- read.delim(yourfile)
5 j2 U8 o/ C1 N( c! `; 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), ] " ] n4 j! K* P" h) I; 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), ]
- x- [- v0 g: I; ]: t dat[i,SST] <- mean(hdf$var)
0 ~, X7 Y" ~" _) p/ T } 3 n1 j3 |% j# K: r% g/ b4 T
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
; ?' r* y7 k0 {/ _- \- N7 j# U 
, V8 ~% ^: _* U! O1 C, I: L' {
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。 
: C6 A+ |# Y- v" F6 z5 ?6 j* o: \) h5 Q
- t) @& f! N6 V; P/ k" X$ f
5 V1 i/ R" a. q0 X) K1 G
$ j# u- a; [! d7 o& M | 
