 使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP）
1 r5 V. o$ M" w+ F7 t Ocean Productivity（http://sites.science.oregonstate.edu/ocean.productivity/index.php）是一个众所周知的海洋生产力数据库，我们经常从中下载相关的遥感数据来用于分析。 $ Z7 m/ i# d2 d0 c! H' m
& U C& [6 S4 z Z7 M7 F
本篇介绍师兄的一个R包，nppr包（https://github.com/chaoxv/nppr）。该包提供了便捷的函数，可以用来下载和处理Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。  ! x5 i# }/ o' T4 }, d7 o5 M0 ~2 u
安装nppr包
. m9 J* z9 r- y 可在github（https://github.com/chaoxv/nppr）获取nppr包。
9 q4 V0 G: F: S' ~$ C& D8 O6 f% a( L #下载nppr包
8 ?, T. U* k. @3 ?' p: F0 H5 m #install.packages(remotes)
" o; N. s5 M2 [' O/ l2 ` remotes::install_github(chaoxv/nppr)
+ S( m# p d% }( l' n$ e$ Y #加载相关R包
$ }) a- r0 o! C1 Y6 Z1 \6 _% B library(nppr) : F3 l1 G5 m0 w6 O+ M
library(RCurl) . n. O3 A, B/ _3 x
library(XML) 8 W% y! S3 W# q8 k/ o
library(R.utils) 6 X4 R% L* `# F# m c I) e
library(tidyverse) % D0 l5 R$ A5 ?$ v8 q5 e$ P: f
library(lubridate)
; j K7 R# A t9 f8 k 使用nppr包下载海洋遥感数据 , `3 |. N* m5 w2 R7 ~9 |+ N3 m
nppr包内常用函数如下所示，get_*等函数可分别用于下载Ocean Productivity的海洋表面温度（SST）、光合有效辐射（PAR）、叶绿素a（Chl a）、净初级生产力（NPP）等遥感数据。 8 {* h' X8 G+ `/ M; W
' W) j% [: W2 u" Y2 h9 ] 以下以获取海洋NPP遥感数据为例作个演示。 8 Z3 u% u. t( d: h, ?4 B
#创建工作路径 8 s8 V6 I9 L, x# }8 G
yourfolder <- F:/R/nppr/vgpm
& D2 R1 Y" E6 s( L( T dir.create(yourfolder)
2 ^7 U3 L' F" D" B' v #以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm ; S+ H; W4 J( T5 [
get_npp_vgpm( 4 y% X' ^# U1 F
file.path = yourfolder, , X- b( q% t+ l" \3 w# A1 b& }9 N
grid.size = low, #指定low或high可更改空间分辨率
& L1 P9 q$ H) M: l: ~/ n l time.span = monthly, #monthly代表月平均，dayly代表8天平均 - U, `; K- c: H6 E: X- W
satellite = MODIS, #选择卫星
* X! n c7 Y' I6 c9 }0 r mindate = 2016-01-15, #指定时间范围以下载遥感
( h( q% Y ]- u. }+ h maxdate = 2016-03-15 + x1 Q8 D1 E$ C# \, d3 K
)
7 l2 d0 \; b( H, u( Q! v 
% U& M) F" @9 m2 _
在这个示例中，我们使用nppr包下载了来自Ocean Productivity（http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php）的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 + y6 a6 t4 N( ]# g: }
使用nppr包进行遥感数据格式转换
/ r3 C3 N- c4 H1 t5 ~ 如上所述，下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法，可将hdf格式转为常见的数据框格式，便于我们后续操作。 # Y$ k& d! t8 D) J ~! X2 A+ b
#将hdf文件转为常见的数据框格式，例如将上述下载的2016年3月的月平均NPP数据做个转换 % N. ~( H$ S2 Q# r7 e
yourfile <- paste0(yourfolder, /201603.hdf)
1 F! l: w9 k( d- } vgpm <- read_hdf(file.path = yourfile)
5 m8 `1 z( s1 M) \7 M2 T5 j9 ^9 D head(vgpm) 0 ]; m8 R' H" ] U) c1 e" D
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
# z8 D% H) J/ ]. ~ 
: u% q, f2 ~+ V: W
转换后的数据框包括三列，分别是经度（lon）、纬度（lat）以及当前日期内该经纬度海区的NPP（var，单位mg C m-2 d-1）。
; @" q6 d0 K+ o. q8 O1 ` 使用nppr包匹配目标经纬度的遥感数据 7 _2 L6 z# D" C
默认情况下，下载的遥感数据是全球海洋的。nppr包同样提供了相关函数，便于我们从中提取特定区域的遥感数据，如下所示。 ; X( b! m! g# b
#获取指定日期和经纬度的遥感，例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP 4 l8 y/ |* l8 A7 w& K! \
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
" Z4 ^: f) h0 P7 ? #或者同时指定多个数据，不再多说 3 |4 o+ h1 Y- G3 ^
mydat <- data.frame(
6 O2 f+ K& p0 C3 ~ lon = c(120, 112, 116),
# K5 B0 B7 K# y: C# D# q lat = c(17, 15, 18), 5 m4 ?6 A/ Y5 p" t K% G. ?# |
date = c(2016-03-04, 2016-03-07, 2016-02-04)
! C4 b& d5 \3 V" u ) % A% k2 y# |" b
match_df(mydat, file.path = yourfolder)  ; y* L' S# L" {
绘制遥感地图
4 ~$ b* m% ~4 ` nppr包的函数geom_oce()可以用来绘制地图，例如我们作图展示来自遥感反演的NPP分布。 ) I) W5 L4 W u& p
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
( ~; _+ e; c- ~ #我们仍以该数据为例作图，展示中国南海2016年3月的月平均NPP
) o( P+ Y9 Z4 k5 G! b# g library(viridis) # v l1 h: b, E; T" W) X7 n
library(ggplot2) g5 x; `+ [3 @
ggplot()+
" j: J, G, V, h9 _+ @4 f geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
& w: n, \1 ^" q% i2 D5 ~ scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
) }# d. U- i! t5 v labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*))) 1 B8 u/ `9 F( j
% i* B& i# U9 H7 X
根据时间和经纬度列表匹配遥感数据的批处理
4 |3 ^8 [) ]1 b$ F5 C7 Z 实际情况中，经常需要对来自不同时间不同经纬度的大量站位匹配遥感，以下提供了一个批处理（不过这是自己先前瞎写的，然后一直偷懒一直用，俺也不知道写的对不对......写在这里仅为方便自己复制粘贴，大家慎用......）
; x8 f' b+ p& R5 Y, f5 N- v 将待匹配的站位的经纬度、日期信息整合在一个文档中，如下所示的这样（本示例命名为“data.txt”）。
* z5 L1 B; ~( Y/ d7 w) t: H 
0 `! ^( A( T m4 f
随后在R中读取该文件，设置一个循环，依次读取日期信息以下载当前日期的遥感（如月平均或8天平均的SST、PAR、Chla或NPP等）。并再根据各站位的经纬度，从中匹配该站位附近的数据（比方说以0.1°为网格进行匹配，并将网格内的数据平均）。
" a! a6 X& X! v5 p) j5 d ##如下以匹配SST数据为例做个演示 + T# R5 {8 `8 f( X. T. W
dat <- read.delim(data.txt) 9 Z9 m7 {3 ?: X/ g
Date <- unique(dat$Date) #获取日期
- F! i/ q$ {( \$ y yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据
o; s5 n1 y! T0 v dir.create(yourfolder)
1 s9 H$ d2 P5 v. G8 s; d" w #通过循环依次获取各日期下的遥感（本示例以下载8天平均SST为例）
3 U6 x4 |! f+ d3 Q for (i in Date) { 4 C9 f9 u+ E# `- E! l. O
yourfolder <- paste(getwd(), SST, i, sep = /) * a4 h) N& _7 l$ j6 J
dir.create(yourfolder)
# Y; Y. g9 j8 ?' D/ W) p get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) # b4 `% q3 p" c8 [6 L
yourfile <- dir(yourfolder)
: P1 w0 ]6 M# e3 ^0 I/ Q hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) # W; F" R2 X8 u8 S% m$ K: U
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)
; x% |' {! T" X+ T! y } : J& D6 X, B+ o8 L F
#再根据列表中各站位的经纬度匹配当前日期的遥感（本示例计算0.1°网格内的平均） . s; o. t/ e/ e6 I
for (i in 1:nrow(dat)) { 0 |! K* A$ u4 S: K
Date <- dat[i,Date] # k/ d% W! M( a" V
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) 6 c* `5 }& {6 e) f/ M
hdf <- read.delim(yourfile) " P' r$ q8 p6 p5 T$ 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), ] * M6 h3 q9 o) C8 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), ] Y7 K4 _- y' i( T! v0 g
dat[i,SST] <- mean(hdf$var) + [8 \9 t. q( Z$ D
}
. {% ?6 b$ Q8 Y2 Q' n0 K3 J) l9 D9 @ write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) / d8 p6 k/ a5 \! _, M
" \$ ^$ o% X" u( _2 C6 H 输出列表的最后一列添加了匹配的遥感数据（本示例匹配了SST）。  # w+ b$ Z d+ ~ z4 X% k1 K+ [! s' Z; K3 J
2 F% `, c( ]( K
6 v; p: `- u0 k0 g
& L' E. D9 }( `& O( i7 S1 e
' _/ Q4 s( U7 V @: E I | 
