使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP)
. _% Q" J2 o* |2 O Ocean Productivity(http://sites.science.oregonstate.edu/ocean.productivity/index.php)是一个众所周知的海洋生产力数据库,我们经常从中下载相关的遥感数据来用于分析。 1 v, s! [4 S+ \6 @

$ i! u3 M# h6 z. r7 A) E* p 本篇介绍师兄的一个R包,nppr包(https://github.com/chaoxv/nppr)。该包提供了便捷的函数,可以用来下载和处理Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。
" z& O/ z% X5 |( O0 v 安装nppr包 ! B$ a G6 S* t# X" @
可在github(https://github.com/chaoxv/nppr)获取nppr包。
2 w4 B0 c- X4 v- B, s+ g: E #下载nppr包
2 t$ X* r& ^7 Z4 C #install.packages(remotes)
6 s4 r# V1 C9 Y& d4 ~4 [+ m remotes::install_github(chaoxv/nppr)
# ~& _* t0 F$ h" { #加载相关R包 . Y* T0 F& r! j# j8 l; P
library(nppr)
! f) o" ]/ F1 h library(RCurl)
6 J: G( N; E0 P library(XML)
) e! \+ O* F" h. z$ Q2 R, X: ]" x library(R.utils)
& Z. W& ?1 K* z! J: P( x library(tidyverse)
" z7 B5 C& C9 m; j/ D: \9 E library(lubridate)
, T' ?. M8 i$ U) Y* } Z, s+ R! d ~ 使用nppr包下载海洋遥感数据 ) j7 f6 G) u8 t) X: B
nppr包内常用函数如下所示,get_*等函数可分别用于下载Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。 6 i' D3 h/ q7 b: x

7 S6 O; q: S+ r" O: h7 [3 w 以下以获取海洋NPP遥感数据为例作个演示。
3 R# D- X7 ]& n/ x$ N #创建工作路径
5 M+ s- O6 o5 U7 k/ I4 M9 e yourfolder <- F:/R/nppr/vgpm * C/ j. T7 H' L
dir.create(yourfolder)
7 [4 h& ~, z: o5 E/ I #以VGPM(NPP的一种遥感算法)为例做个演示,详情?get_npp_vgpm
" N' I( }& b. R& G- r4 r get_npp_vgpm( * D' f# ]' [/ w' p/ q
file.path = yourfolder,
- n% E' z) `' n7 y2 K% c grid.size = low, #指定low或high可更改空间分辨率
, |7 N7 ~- t0 d4 x p5 }$ Z _3 p time.span = monthly, #monthly代表月平均,dayly代表8天平均 * x' q/ T' N5 \3 X1 J" Y
satellite = MODIS, #选择卫星
( g: ~4 D/ q, Z1 L. y& }: Q mindate = 2016-01-15, #指定时间范围以下载遥感
5 x4 W: J+ w) C# c3 N5 e9 O maxdate = 2016-03-15
8 M0 S$ L& Z1 H$ X+ ` )
5 O7 U$ d; G# C: \% r( ^, w  : v4 X: H; A R5 @, b5 z
在这个示例中,我们使用nppr包下载了来自Ocean Productivity(http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php)的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 . _+ {% F0 S* o6 Y3 A
使用nppr包进行遥感数据格式转换 $ c& S9 ?# L1 t8 ]( e) P! k
如上所述,下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法,可将hdf格式转为常见的数据框格式,便于我们后续操作。
( x' ]5 ]! _4 P7 L #将hdf文件转为常见的数据框格式,例如将上述下载的2016年3月的月平均NPP数据做个转换 + _- p+ z$ E6 d# \
yourfile <- paste0(yourfolder, /201603.hdf) $ o: s% y8 B9 d" W
vgpm <- read_hdf(file.path = yourfile)
) H! E1 O3 t% _- I head(vgpm)
/ {8 d% J" k. u+ z0 Y write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
- @% z+ |6 \6 [  / W3 o" ]: k- I5 s* k; J
转换后的数据框包括三列,分别是经度(lon)、纬度(lat)以及当前日期内该经纬度海区的NPP(var,单位mg C m-2 d-1)。 2 Y$ v& U" O+ n \8 ]- R5 a
使用nppr包匹配目标经纬度的遥感数据 : h/ X6 b2 V. V8 A
默认情况下,下载的遥感数据是全球海洋的。nppr包同样提供了相关函数,便于我们从中提取特定区域的遥感数据,如下所示。 ; n) R0 X6 X8 l' r3 s" k& [ i6 e
#获取指定日期和经纬度的遥感,例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP / m. H* U) M8 h( A& J- r! z1 I
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
3 `0 h% a7 g$ N #或者同时指定多个数据,不再多说 & f1 M; d) ~/ G, s k: I0 Y+ ~
mydat <- data.frame(
8 r( ?! n8 Z! _- u- b9 n; u lon = c(120, 112, 116),
8 p3 A+ X9 P* e; C! ` lat = c(17, 15, 18), : _ ]! M: @; j2 |% e& V C0 S* C
date = c(2016-03-04, 2016-03-07, 2016-02-04)
9 T" w( y7 E7 j& i9 B2 L ) ( p4 _( N0 ~' w/ l4 C9 \
match_df(mydat, file.path = yourfolder) ( G! @& D. @& Q/ S
绘制遥感地图
. o* l0 `. ^$ [% U nppr包的函数geom_oce()可以用来绘制地图,例如我们作图展示来自遥感反演的NPP分布。
3 M# p4 v$ ?1 n# }/ [' | #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式 1 [, U9 G. }; H4 n3 l3 G
#我们仍以该数据为例作图,展示中国南海2016年3月的月平均NPP 9 F4 t4 z0 R9 S- }# b* n C3 g
library(viridis)
4 F/ l9 M4 Z0 a: W library(ggplot2) 4 w1 y* F# l2 l1 C' E
ggplot()+
' Q% z/ L4 n% K4 C2 q3 R) ] geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+
/ i. g6 h: S6 I( c scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+ / J! a7 A5 {" j( G" i' \
labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
* X* E, E0 y+ C+ H' \  9 \% h8 O$ B3 p0 Q. Y4 e$ g8 h
根据时间和经纬度列表匹配遥感数据的批处理 / F& [ g( B T2 C1 P) i J
实际情况中,经常需要对来自不同时间不同经纬度的大量站位匹配遥感,以下提供了一个批处理(不过这是自己先前瞎写的,然后一直偷懒一直用,俺也不知道写的对不对......写在这里仅为方便自己复制粘贴,大家慎用......)
9 j2 W1 o' N' _ U1 c5 ]+ T 将待匹配的站位的经纬度、日期信息整合在一个文档中,如下所示的这样(本示例命名为“data.txt”)。 7 x& X' b) [' l# X, D) Q

) `/ J* c- t7 E- b 随后在R中读取该文件,设置一个循环,依次读取日期信息以下载当前日期的遥感(如月平均或8天平均的SST、PAR、Chla或NPP等)。并再根据各站位的经纬度,从中匹配该站位附近的数据(比方说以0.1°为网格进行匹配,并将网格内的数据平均)。 ! I) U& |/ {3 F* U0 W/ P
##如下以匹配SST数据为例做个演示 0 P, j1 K6 N; ^( }" g9 X
dat <- read.delim(data.txt)
" D0 b. }: H _, U% O5 C9 V' { Date <- unique(dat$Date) #获取日期 3 N, V2 Q1 ^6 M/ B( G
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 5 W5 C" N/ ~$ E3 X" N8 s3 ]6 K& \7 o! P
dir.create(yourfolder)
1 p2 ?% Z2 x5 F #通过循环依次获取各日期下的遥感(本示例以下载8天平均SST为例)
9 V3 h9 m6 ]7 L# Y4 I for (i in Date) {
. c1 t* h' b" {# l0 b7 v yourfolder <- paste(getwd(), SST, i, sep = /) * `8 q2 o4 C) V0 m1 F; n% }+ S
dir.create(yourfolder) ( O9 j2 f: o: n. P' |% e
get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) ' V u$ d1 f) P7 K' A* P. f" E+ k
yourfile <- dir(yourfolder)
, R& E$ k% y! R, j hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) 2 g8 ], n% x n3 A& E! D Q% q
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) , A7 ?/ J* y0 u5 C- s O! q
} " d5 a3 {3 c* h' ?% c% [
#再根据列表中各站位的经纬度匹配当前日期的遥感(本示例计算0.1°网格内的平均) ) d3 o- q" O5 [5 c* j3 B
for (i in 1:nrow(dat)) { * n" }) x8 @) ?+ \
Date <- dat[i,Date] 3 N6 M2 U9 j5 g# j7 {/ @
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )
( W- p. c. J# @' I5 T3 z( A hdf <- read.delim(yourfile) 6 i0 `% N) M4 t" V
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), ]
! H, _1 `5 Z0 {: M$ t' y+ E$ A" 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), ] 7 y5 U$ J/ q% Z$ I/ W# m$ k r% o
dat[i,SST] <- mean(hdf$var) 4 j0 Z% K$ t, E- E. X
}
% d% p( K6 b8 _ write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) % j$ d% T3 J% p* t; _' F' J
 7 B" [$ T( Y$ d
输出列表的最后一列添加了匹配的遥感数据(本示例匹配了SST)。 ! n; [: r, c Y# @. c& k
! d. t- ?/ b$ {) Y5 J
/ L6 m) H, ^ K3 M5 `/ G, o5 D1 w) u! {. a% R: y" `& s( Z! v1 C
7 _9 @1 e- ~: R6 v& `& H" r |