 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） 5 \7 I# O9 H9 Z5 ~0 K8 d$ X8 Z
Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。
2 w) o- {+ \: W* T% j% e" M 
7 l' m1 v, e8 ` [2 W- I; A2 _ 本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 
G0 r( ^& D& q7 _7 k1 j 安装nppr包
& q% f7 f* C" q! x& Q {4 L, V6 g8 x 可在github（https://github.com/chaoxv/nppr）获取nppr包。
; t& G) `7 p) O: ?. x* h" ~& O/ W #下载nppr包 # f1 j% h6 g1 n* c4 ?$ y. `
#install.packages(remotes) - N' _- x, R4 {4 |; Q
remotes::install_github(chaoxv/nppr) ) N( E: I. ?1 {
#加载相关R包
+ A/ M" }. s; l) \ library(nppr) 8 Q V1 M9 m4 s1 ?% G( a+ P
library(RCurl)
) [0 t- g5 R! G% j library(XML) / Q& E! \. [ g1 Y; w5 S% a
library(R.utils)
0 z/ j. `# L0 u library(tidyverse) - d4 |* g( d0 k! [& k5 r, ^# _
library(lubridate) * s. M: p( {3 ~2 T) B0 ~/ K! d6 K
使用nppr包下载海洋遥感数据
+ N$ d- ]: z# h nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 % c2 X" S6 }: {9 W9 H
! U; c0 h; T0 x i* c( s8 ?" M 以下以获取海洋NPP遥感数据为例作个演示。 0 g' n" J$ _. f( _: U1 [
#创建工作路径 & f: G- ~& f! c
yourfolder <- F:/R/nppr/vgpm ( w" p& K' b( t8 o
dir.create(yourfolder) $ l: \2 Z0 F+ `) _" E+ k& k
#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm
* Z k+ n3 ], c" U' g7 c" z get_npp_vgpm(
. Y+ @; `. I9 Q q file.path = yourfolder, + C1 s3 V. t+ g5 S F' B1 o
grid.size = low, #指定low或high可更改空间分辨率 ) V# r. e7 M0 e% O" A
time.span = monthly, #monthly代表月平均，dayly代表8天平均 ; n1 `! y4 t- R+ V7 g9 @! h
satellite = MODIS, #选择卫星
6 i; j& l0 Y& p8 X" S mindate = 2016-01-15, #指定时间范围以下载遥感 3 {" q: [5 A) v% a7 g9 f
maxdate = 2016-03-15
# W p- ^! i$ p7 ]- g ) ) z. h% z8 m& p6 Z' r) \
3 i3 A3 Y* O5 J
在这个示例中，我们使用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 A C! Y. h* X2 H! F2 V 使用nppr包进行遥感数据格式转换
# i6 ^- c1 G- G' j9 z 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。
, ], J4 } L( H' z- c& g2 f( f #将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换
6 C8 K9 Z- c/ `4 X( c/ P yourfile <- paste0(yourfolder, /201603.hdf) 8 E g) b+ ^ f' K
vgpm <- read_hdf(file.path = yourfile) , w8 ~; q" z; U0 m% u- e6 I
head(vgpm)
! q [0 a6 Z7 \) z, F, i8 }" X4 w a write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
+ ^6 G3 X$ N: I 
3 p3 e) ?7 N# X/ c
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
, m: r9 A% @- H; l7 n4 w 使用nppr包匹配目标经纬度的遥感数据
$ }1 d/ r) ^; ~; _2 U; I2 q6 z* S0 j 默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 $ h. Q' G |8 P. v) g
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
/ f, ?! q5 t; A+ j/ Y match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01) 4 R; I7 L1 M$ D( F
#或者同时指定多个数据，不再多说 Q9 @; Q& r# X% d [' d' S
mydat <- data.frame(
' v. M+ N. n3 l lon = c(120, 112, 116), 7 A; w+ E* A( m- ?4 |" B
lat = c(17, 15, 18), 7 d9 N& M) e' p. q. T4 C( ^% A2 d0 _
date = c(2016-03-04, 2016-03-07, 2016-02-04) ; g0 N2 |" E1 Q6 j7 H5 f7 Q
) ! I$ e! y- G& ?. Z/ Q; W
match_df(mydat, file.path = yourfolder) 
e1 X( m! ]$ D/ o6 c. p# b 绘制遥感地图 % z8 h( j/ @/ U- y" e
nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。
5 e: w- X7 `. v #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 7 q, W) r" s# ~( L
#我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP 8 X9 l7 p g3 ]4 E3 ]7 [
library(viridis)
% y! l/ K& @' M/ K X% E2 ? library(ggplot2)
2 n: z. O" P* ]/ z+ n [: P8 C ggplot()+ 5 ~7 z$ g7 @6 f* n) j4 J
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
0 z% X. Z& ^' ?. p. I scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ 5 J* h( N( P0 w" ^+ h
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) 6 v+ I1 g9 [2 o: [! n1 s' ?
+ M) ^, G ~$ q' t: ~. V
根据时间和经纬度列表匹配遥感数据的批处理 : t8 o! t; ^- C1 Q5 a: X/ y
实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......） 3 t" | O7 g7 o6 g. e
将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。 ( C Z# |" h9 b M' J
. u& r8 M" I. R3 w 随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。 # |9 X; D: ^$ [
##如下以匹配SST数据为例做个演示 + ^$ d; x) m4 W( \
dat <- read.delim(data.txt)
9 F% M' G0 {+ |# _7 s7 Q6 o$ s7 } Date <- unique(dat$Date) #获取日期
8 }' f1 _* ]! n/ h0 y, g yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据
7 j! ]) ]5 i* P dir.create(yourfolder) " K5 a, Z$ |' u* ^$ u! I! v
#通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例） 5 X5 ^$ z. Z* n% q
for (i in Date) {
3 h" R+ i+ j' w: P6 s! [% l( l7 ? yourfolder <- paste(getwd(), SST, i, sep = /) " k. u' r, w8 u& @! j
dir.create(yourfolder) ; X# u! e* `1 ?1 F
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) 2 h9 g8 D7 ^4 {
yourfile <- dir(yourfolder)
. `( f1 e; C8 D+ e; d hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))
4 {, A( n) |5 K( M! \1 x& M write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) 5 Z% i7 E' C- V c( C) I
} 7 Q# F. k9 ^% p1 i, P
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） + {5 K0 Y, ^4 R
for (i in 1:nrow(dat)) {
' F. {& R1 b' d. }3 o2 v( |% n: ^$ Y5 ? Date <- dat[i,Date] + i9 z; p- ]& Z2 N4 A) @( }4 I
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) : ?/ u- B e% i
hdf <- read.delim(yourfile)
, k, E9 O$ u, _& d3 s- _1 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), ] , l) l! O9 s% u6 f% [
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), ]
7 m$ g/ h0 T4 ]5 z: k& z) |* q dat[i,SST] <- mean(hdf$var) 0 x7 N, ^" I5 V+ X9 e y
} : i0 Z' K3 I$ R) M1 W. R9 q" e
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) + W9 T% W+ j. Y9 ^; N
) C9 I+ C9 y2 ~1 r. O( x
输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。  4 Q! J; W; b* t; L4 P7 S& N
- q0 A9 Q9 i( w1 \) n2 {4 F
1 N( W6 W9 O! I) a- [6 o
) Q+ o e+ E5 k
% P: q2 r3 y, [6 k& S | 
