使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP) / V; Y0 {* ^4 E6 p
Ocean Productivity(http://sites.science.oregonstate.edu/ocean.productivity/index.php)是一个众所周知的海洋生产力数据库,我们经常从中下载相关的遥感数据来用于分析。
9 _3 J' k2 ~8 o  ' w2 C5 M2 }+ n. V
本篇介绍师兄的一个R包,nppr包(https://github.com/chaoxv/nppr)。该包提供了便捷的函数,可以用来下载和处理Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。 , l& B/ `* r" A1 L4 L
安装nppr包 ( W. X: [' ]% o! r5 z# U
可在github(https://github.com/chaoxv/nppr)获取nppr包。 # j- a. U* C) |; E, D9 O: O; d
#下载nppr包
7 @$ _) Q# W! u. }8 V #install.packages(remotes)
" J& n- F6 \8 l7 G% E. v remotes::install_github(chaoxv/nppr) $ \% t1 ~1 ~* }# _4 k* Y
#加载相关R包
) b# P3 }" Y$ r g library(nppr) Q6 H: E! G! l2 b4 y( x+ p" g
library(RCurl)
( |& p# B4 V) g! a& T library(XML)
9 N0 o4 f4 q* s- B$ [& y3 ` library(R.utils)
- r" h" O4 q8 I library(tidyverse) ' v, g* w- M8 m( R: C# H( a
library(lubridate)
+ @1 c0 X% X: T 使用nppr包下载海洋遥感数据
! ?8 c5 W) `( v nppr包内常用函数如下所示,get_*等函数可分别用于下载Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。 1 Z. _ R; j1 B+ y( j# d+ r

4 }4 I) q2 ~" F5 b 以下以获取海洋NPP遥感数据为例作个演示。
x) \, f3 {8 m5 A/ j& } #创建工作路径 . s j! D4 b% Y& f+ S
yourfolder <- F:/R/nppr/vgpm j) W( s3 G8 q1 W$ ]
dir.create(yourfolder) # o( _7 l* H! [8 ^7 W+ ~. q
#以VGPM(NPP的一种遥感算法)为例做个演示,详情?get_npp_vgpm
. u, [3 S2 e; r, A5 O1 E8 t get_npp_vgpm(
/ z" p) z. R4 N% \# }7 v file.path = yourfolder,
5 G1 T' r! r1 O( Z grid.size = low, #指定low或high可更改空间分辨率
9 l* O. T* ?0 n4 X; K0 t/ T( V time.span = monthly, #monthly代表月平均,dayly代表8天平均
# _/ j5 L9 g$ H } b0 b! p satellite = MODIS, #选择卫星
! u: q$ y. [! ^9 x9 p mindate = 2016-01-15, #指定时间范围以下载遥感
8 ]7 {' p8 Q& y/ ]1 O+ S/ z7 P2 ] maxdate = 2016-03-15 / R9 U$ D0 F- ?' B
)
4 a. _) ~, z3 W' p$ U; V% B! f+ Z 
: F1 g& `2 w9 S8 V 在这个示例中,我们使用nppr包下载了来自Ocean Productivity(http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php)的基于VGPM算法反演的全球海洋2016年1月至3月的月平均NPP数据。 . i' ^5 W* e( H* m F, L/ }$ a
使用nppr包进行遥感数据格式转换
! P8 m6 F* V" i3 z! [ 如上所述,下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法,可将hdf格式转为常见的数据框格式,便于我们后续操作。 ?! h' z9 ^' D2 ?" \
#将hdf文件转为常见的数据框格式,例如将上述下载的2016年3月的月平均NPP数据做个转换
( s, ?9 Y! Z/ } yourfile <- paste0(yourfolder, /201603.hdf) 8 s( e; U: o# S1 J6 U& @7 R0 ], m4 |
vgpm <- read_hdf(file.path = yourfile)
2 v, |% f& @: I# p# S8 @6 B5 d/ l8 p head(vgpm)
$ w' c7 w7 j/ p4 Y6 ^ write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)
W" F* w7 d4 m  5 W; ]1 ]; L% D6 A4 G
转换后的数据框包括三列,分别是经度(lon)、纬度(lat)以及当前日期内该经纬度海区的NPP(var,单位mg C m-2 d-1)。
, m% ~/ O" g/ T) H' ] 使用nppr包匹配目标经纬度的遥感数据 0 G. l. W* m6 ~9 L" R
默认情况下,下载的遥感数据是全球海洋的。nppr包同样提供了相关函数,便于我们从中提取特定区域的遥感数据,如下所示。 ! x5 h- w1 D' r
#获取指定日期和经纬度的遥感,例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP
+ d( u: X# \9 n match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
/ k: P- m4 t5 O+ D# D& _4 ] #或者同时指定多个数据,不再多说
4 S( n G; i1 B. v: j mydat <- data.frame(
5 X' l" F. Y7 k k# U ^, ]) ? lon = c(120, 112, 116),
+ F- ~, c9 J% r. n' k5 K; e lat = c(17, 15, 18), 4 C1 Z! {7 k( j# w- [+ s4 M
date = c(2016-03-04, 2016-03-07, 2016-02-04)
5 A4 ^2 B7 W1 I. H ) ) c4 [" s+ h6 E( z. f2 C
match_df(mydat, file.path = yourfolder) ( E! C: p$ h1 d: x& ]
绘制遥感地图
6 }# D" E4 ~6 ?0 K nppr包的函数geom_oce()可以用来绘制地图,例如我们作图展示来自遥感反演的NPP分布。
# g1 ^3 M( F L. {0 K #上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
# z" M$ y% X3 @ #我们仍以该数据为例作图,展示中国南海2016年3月的月平均NPP
+ g' }+ \' N, Y6 W. L library(viridis) & @9 `; {; I$ X- H) T. Y, W
library(ggplot2) 0 D7 t2 H8 r8 Q3 m% @3 P; y
ggplot()+ ( W3 O: ?7 `& j8 Z
geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ 0 A3 O* j: Q' d- l
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
! m0 f7 P' ?$ [' @+ m5 v labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
& I: l6 {: Y* z$ q9 w  " O2 W6 [6 b/ ^/ L& |" ?$ Q* t
根据时间和经纬度列表匹配遥感数据的批处理 5 e7 y( s& m Q6 b
实际情况中,经常需要对来自不同时间不同经纬度的大量站位匹配遥感,以下提供了一个批处理(不过这是自己先前瞎写的,然后一直偷懒一直用,俺也不知道写的对不对......写在这里仅为方便自己复制粘贴,大家慎用......) # ~. \3 k' S/ @+ h
将待匹配的站位的经纬度、日期信息整合在一个文档中,如下所示的这样(本示例命名为“data.txt”)。
8 @+ t$ l8 [& A! o* } 
1 A* j2 h. f7 ^7 C# Q) _3 N2 M1 u2 h; u 随后在R中读取该文件,设置一个循环,依次读取日期信息以下载当前日期的遥感(如月平均或8天平均的SST、PAR、Chla或NPP等)。并再根据各站位的经纬度,从中匹配该站位附近的数据(比方说以0.1°为网格进行匹配,并将网格内的数据平均)。
; k+ I; B8 V/ r2 ` ##如下以匹配SST数据为例做个演示
+ N8 ~# _' U5 P dat <- read.delim(data.txt)
; q# q% Y. h. f6 N& j/ v Date <- unique(dat$Date) #获取日期 " {2 q' I* w3 c4 O; V+ G2 `
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 + ]4 C" U2 }1 J, ]9 G1 g4 ^
dir.create(yourfolder)
! N& Z/ y% t4 B% J! C2 N. j #通过循环依次获取各日期下的遥感(本示例以下载8天平均SST为例) % |8 b U* l2 e0 W
for (i in Date) { ) q4 `: v5 I) x
yourfolder <- paste(getwd(), SST, i, sep = /)
; r5 @3 s9 U9 C% j. @, u$ I" a dir.create(yourfolder)
4 w/ ~ M% {/ m, \0 ]: |' M- l get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i) 1 T) @& G) q, _2 `9 r+ ]8 {4 W8 B. M6 x
yourfile <- dir(yourfolder)
* p2 G) Y& t9 E) i; W9 ]' t$ W0 m hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) " Y, D/ A: R* {* Q1 e# B4 A
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) & e9 X8 ~' e$ N2 b* C7 j, n' C a
}
, O, C' H( \! \8 q5 Z% V( W, k #再根据列表中各站位的经纬度匹配当前日期的遥感(本示例计算0.1°网格内的平均) $ o# p7 z1 ~$ C0 ? t1 w C
for (i in 1:nrow(dat)) { * t* I- z4 N7 \& S( F( u2 s
Date <- dat[i,Date] ( t0 y, T+ S1 z0 y
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )
& `' |. _9 D6 e/ R4 G% m, x B hdf <- read.delim(yourfile) ' ~: n6 K- V& b* |
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), ]
# D2 u# R- a) l 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 o9 F. j' c# f4 g) o
dat[i,SST] <- mean(hdf$var) - |' m" T; O" H; j2 K+ o) R
} " n1 A0 D; _) S) Q6 |7 e
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)
. Q! j* k- R& P5 B! k+ ^ j3 i  % m( p. s. o( S1 i
输出列表的最后一列添加了匹配的遥感数据(本示例匹配了SST)。 2 H- O+ M2 X% i p& e2 L
. V" S8 u l2 l7 Y3 K( F* u: P5 p. d' l
# z9 o6 O) O( A
" l6 e5 P5 O. r. q. I/ f4 u |