使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP)& W! }/ x j/ P/ g0 S8 Q. f
Ocean Productivity(http://sites.science.oregonstate.edu/ocean.productivity/index.php)是一个众所周知的海洋生产力数据库,我们经常从中下载相关的遥感数据来用于分析。
; _/ C4 H& R- n5 o! @9 d
1 Z4 r; b) Q$ V) h4 X Y 本篇介绍师兄的一个R包,nppr包(https://github.com/chaoxv/nppr)。该包提供了便捷的函数,可以用来下载和处理Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。 * x' k/ j. ]) A5 { D
安装nppr包
0 _- [ s: U3 X9 V l& M 可在github(https://github.com/chaoxv/nppr)获取nppr包。 9 M4 I& `* O0 ^: ~$ c/ b. P
#下载nppr包 ; y- i/ y& Q8 \/ }1 \( Z
#install.packages(remotes)
1 X$ b9 l' q( ] Z1 x6 t! V remotes::install_github(chaoxv/nppr) . |; S4 v/ L" \3 G1 E; e
#加载相关R包
4 {* V0 z# O1 [+ d* L4 f R3 v library(nppr)
4 B: |. `1 w- a' d! z library(RCurl) . C) q; J6 U3 s; c
library(XML)
; K5 z+ O& i. z library(R.utils)
7 ]( E6 K4 y/ r+ z" C! s7 j library(tidyverse) + v$ ~2 |( r: O+ ]; v( O7 n
library(lubridate) - U2 M) p; K5 J: Z4 e1 m9 V
使用nppr包下载海洋遥感数据
2 G1 g9 `8 h' O6 _ nppr包内常用函数如下所示,get_*等函数可分别用于下载Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素a(Chl a)、净初级生产力(NPP)等遥感数据。 , ^: [6 c5 w6 l, s
9 w$ y" ~; n/ ^. O1 R5 \% o1 n! k 以下以获取海洋NPP遥感数据为例作个演示。
$ p" b$ O/ ?9 ~1 j #创建工作路径 # z3 v# c) a6 z$ R+ j' }3 O1 U: I
yourfolder <- F:/R/nppr/vgpm
6 q O; ~: ~3 u* W7 n7 c dir.create(yourfolder)
' B" b5 {* {4 R2 f+ X0 ^, E! d/ W #以VGPM(NPP的一种遥感算法)为例做个演示,详情?get_npp_vgpm _- N9 W4 s/ N! }
get_npp_vgpm(
+ j R5 `. a8 P! @4 B: ^3 ? file.path = yourfolder, 4 ~4 \5 _% \. t! X8 ?
grid.size = low, #指定low或high可更改空间分辨率
5 T. ?8 F6 z5 j$ N3 p/ D, z time.span = monthly, #monthly代表月平均,dayly代表8天平均 ( n0 F* J# `: Z- @) J
satellite = MODIS, #选择卫星 : Y& v" ?6 C( B/ Z! E# }- j
mindate = 2016-01-15, #指定时间范围以下载遥感
% C6 t- v L7 V+ N) m0 g% ]4 X$ j' U maxdate = 2016-03-15
' h% ~2 v2 L ^0 f6 b- I1 }. k) O ) ( U7 g: l- z+ @7 o# f
( C. |8 l! A8 y, @3 }# y& 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数据。 2 o2 O) ?% ?6 `8 g& y
使用nppr包进行遥感数据格式转换 / @/ q! f0 d/ q3 B- {% e* q! s- h
如上所述,下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法,可将hdf格式转为常见的数据框格式,便于我们后续操作。
4 L% f5 ]; N* l! T y8 C. t4 _ #将hdf文件转为常见的数据框格式,例如将上述下载的2016年3月的月平均NPP数据做个转换
9 @3 D: q, Z$ Z- U7 D yourfile <- paste0(yourfolder, /201603.hdf) 9 e1 \1 R. d. l3 R2 Q/ ^7 H3 }
vgpm <- read_hdf(file.path = yourfile)
/ W5 d4 N- Z9 G( K head(vgpm) : i4 ^. o2 ]5 y5 Q" g
write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE) ( x9 O7 m# ]3 e9 }; `9 X& S
: k: r6 w) s) I5 g! Y
转换后的数据框包括三列,分别是经度(lon)、纬度(lat)以及当前日期内该经纬度海区的NPP(var,单位mg C m-2 d-1)。
/ T+ N& ?4 j: Y Q) L) \; U 使用nppr包匹配目标经纬度的遥感数据 ' u& D4 E- N4 |' V$ e
默认情况下,下载的遥感数据是全球海洋的。nppr包同样提供了相关函数,便于我们从中提取特定区域的遥感数据,如下所示。 1 q( S# R+ P3 H0 u7 d* [' }0 s
#获取指定日期和经纬度的遥感,例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP & J6 X! R3 ]# Q1 d" K1 m
match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)
+ ?) P+ x' O) w7 s7 X) Y- _ #或者同时指定多个数据,不再多说 ) D; y2 y4 w } o* y
mydat <- data.frame( 4 C! E1 P' \6 i: u$ n+ [. L6 M
lon = c(120, 112, 116),
+ t$ f4 E* p7 Y4 T5 v1 @0 @0 k' K lat = c(17, 15, 18),
) G" [ Q- i) B# J5 d) E# R1 n& @ date = c(2016-03-04, 2016-03-07, 2016-02-04)
3 B. G( B, a. J% M )
1 C5 z5 b, ~* b0 ^3 T4 p match_df(mydat, file.path = yourfolder) * F/ w, g+ ?5 k, p" Y4 N
绘制遥感地图
j8 M/ v2 P+ h6 p) G+ p9 L nppr包的函数geom_oce()可以用来绘制地图,例如我们作图展示来自遥感反演的NPP分布。 ! l# U" ~# m1 \: O- a$ v" b
#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式
3 l7 T% u# C; |$ q3 O #我们仍以该数据为例作图,展示中国南海2016年3月的月平均NPP
: B1 `, U2 Q, b% T6 K# c library(viridis)
# l, h- R2 x) z7 d; @1 q# u/ i library(ggplot2)
, l4 u" T7 k$ B& P0 l1 X ggplot()+
; J, d! Z* f1 D) ? geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+ / C0 p* J" G. D. H9 ~. ^/ j
scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+
9 y' ], [* e5 [- g( Q' K labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))
& ^9 d8 u- z* H
4 m2 a& y f6 L! L& ^+ d: V3 P8 j 根据时间和经纬度列表匹配遥感数据的批处理
8 C/ X9 a6 S# ]4 f/ u0 k' }* t 实际情况中,经常需要对来自不同时间不同经纬度的大量站位匹配遥感,以下提供了一个批处理(不过这是自己先前瞎写的,然后一直偷懒一直用,俺也不知道写的对不对......写在这里仅为方便自己复制粘贴,大家慎用......)
1 T2 f& g+ P2 j# K 将待匹配的站位的经纬度、日期信息整合在一个文档中,如下所示的这样(本示例命名为“data.txt”)。
7 T4 _7 R/ i( p8 [7 T 4 Y* V8 p* Y( w: t; K
随后在R中读取该文件,设置一个循环,依次读取日期信息以下载当前日期的遥感(如月平均或8天平均的SST、PAR、Chla或NPP等)。并再根据各站位的经纬度,从中匹配该站位附近的数据(比方说以0.1°为网格进行匹配,并将网格内的数据平均)。 * M7 h; z& N1 W2 _/ A$ K
##如下以匹配SST数据为例做个演示 % d, h+ h+ `, J% @
dat <- read.delim(data.txt)
7 s7 t& `% E2 c* G- z. F; e* N Date <- unique(dat$Date) #获取日期 9 c' ?* I) l4 q+ ?) U8 z# t- _$ m
yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据 2 b; f, R( B Z% K, e
dir.create(yourfolder)
b6 r% B0 C+ P #通过循环依次获取各日期下的遥感(本示例以下载8天平均SST为例)
- T7 [ e: K) I$ S/ C for (i in Date) { 8 ]. x! Y) ^# O O' u
yourfolder <- paste(getwd(), SST, i, sep = /) 1 v6 R* c( d# D. j
dir.create(yourfolder)
" R+ r- x/ u& O7 u# b" y get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)
7 F& H. r5 f% j4 a yourfile <- dir(yourfolder)
9 D8 `& @9 L6 o" g) B3 A/ s2 r hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /)) 5 X2 g P3 I. V L
write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE) * @5 t) E, T- h# Y5 x# q3 p
}
- x- y. k0 Z/ T6 O6 ~3 E+ O( g #再根据列表中各站位的经纬度匹配当前日期的遥感(本示例计算0.1°网格内的平均) ; b5 E& d: _: Q! J! X" M
for (i in 1:nrow(dat)) {
o" _: f2 q# {2 l. t8 n' {/ @ Date <- dat[i,Date] : u5 n9 t- N" v2 T
yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = ) 5 a: z' x6 s5 \; E
hdf <- read.delim(yourfile) 8 W' q# i# g6 T" l2 u
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), ]
. }6 g! T0 ~: o2 d& t0 w9 z 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), ] y; X- _5 V) X7 V4 C. k
dat[i,SST] <- mean(hdf$var) 5 ~- ]2 I7 x$ n1 h
} 2 _ d8 z0 a% ~7 {" {1 F: |) L
write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE) + e% L; z! E7 p' ?4 \$ v9 V
9 ?& L& s5 m$ ?
输出列表的最后一列添加了匹配的遥感数据(本示例匹配了SST)。 , o% K/ F; S2 q7 G4 h
2 w) Z# T' V/ c3 C
+ F( K1 f S) k* H, _. V
) q9 e- p9 j4 s1 V& E
" E4 V+ o3 n) i8 D7 z3 e
|