使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP) - 海洋遥感数据处理

[复制链接]
使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP) 1 J% k/ [# X! w

Ocean Productivityhttp://sites.science.oregonstate.edu/ocean.productivity/index.php)是一个众所周知的海洋生产力数据库,我们经常从中下载相关的遥感数据来用于分析。

9 a7 e3 h s. a3 }

- E+ X0 E' q3 m+ R& o2 X

本篇介绍师兄的一个R包,nppr包(https://github.com/chaoxv/nppr)。该包提供了便捷的函数,可以用来下载和处理Ocean Productivity的海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素aChl a)、净初级生产力(NPP)等遥感数据。

* O0 V( R; U6 v7 d$ Q- k2 Y

安装nppr包

; l* O4 ~8 w, a! ~) I1 u! I9 O$ \

可在githubhttps://github.com/chaoxv/nppr)获取nppr包。

- f5 G6 i/ ~/ p0 a: P( D

#下载nppr包

5 G4 a* W2 w% H1 O6 b! n9 b

#install.packages(remotes)

+ G: P- W+ Q2 u# r4 |, `/ {$ k

remotes::install_github(chaoxv/nppr)

7 v" ?5 f4 g6 C+ g; k

#加载相关R包

$ _2 x2 x/ m+ N( W+ I$ _% k

library(nppr)

?- z; I p" `! M

library(RCurl)

! a0 f! V' U, G' U: P6 o( A

library(XML)

, _" n3 ~2 S% v3 O- l9 R$ C

library(R.utils)

W8 s3 Y8 U8 v- s' s) W

library(tidyverse)

: Y2 `4 f( T7 g# Q0 R/ h" j

library(lubridate)

1 J+ u1 ?; q' q) Y. t- V9 [

使用nppr包下载海洋遥感数据

! P3 N# P$ u- w: F

nppr包内常用函数如下所示,get_*等函数可分别用于下载Ocean Productivity海洋表面温度(SST)、光合有效辐射(PAR)、叶绿素aChl a)、净初级生产力(NPP)等遥感数据。

2 t, t4 R. G' {- {- t6 @

- l. {: s) A. a2 S

以下以获取海洋NPP遥感数据为例作个演示。

a, x; ^$ g+ n2 W5 o& u7 v9 ~ c

#创建工作路径

9 }# }* ~& L) f

yourfolder <- F:/R/nppr/vgpm

8 g6 V5 w ?8 D5 }

dir.create(yourfolder)

! S6 ^4 o2 P5 |. U k7 h

#以VGPM(NPP的一种遥感算法)为例做个演示,详情?get_npp_vgpm

5 Z* I% }$ I: T l4 ?5 A

get_npp_vgpm(

8 p G- k! d) C* z2 N! S

file.path = yourfolder,

) U$ D) S8 w: \$ F3 X) V, x6 j

grid.size = low, #指定low或high可更改空间分辨率

/ E9 I; {; n: |7 i8 J+ P

time.span = monthly, #monthly代表月平均,dayly代表8天平均

- i! K4 C2 a" d I7 A7 _

satellite = MODIS, #选择卫星

) N# a; J3 p0 O+ C1 M( z

mindate = 2016-01-15, #指定时间范围以下载遥感

v3 h& o. }0 H# b, E( o) O) A6 h

maxdate = 2016-03-15

# l% T3 c3 h- z6 `/ R- l3 \

)

- v4 l' p7 V! p6 ^) o

0 O$ p1 p" h; U- N

在这个示例中,我们使用nppr包下载了来自Ocean Productivityhttp://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.vgpm.m.chl.m.sst.php)的基于VGPM算法反演的全球海洋20161月至3月的月平均NPP数据。

+ W: ]' }7 |1 j1 \ D x8 a

使用nppr包进行遥感数据格式转换

. A5 k2 H# a7 Z' a) l& P1 X

如上所述,下载后的遥感数据以hdf格式存储。nppr包提供了便捷的方法,可将hdf格式转为常见的数据框格式,便于我们后续操作。

& d. X6 @1 q1 A! I& \! T

#将hdf文件转为常见的数据框格式,例如将上述下载的2016年3月的月平均NPP数据做个转换

( ]0 F6 H# e+ c9 Q- u

yourfile <- paste0(yourfolder, /201603.hdf)

# y" |& [8 C6 t0 D

vgpm <- read_hdf(file.path = yourfile)

2 P2 e. _; N& e) E* s# E+ E+ c

head(vgpm)

9 N1 }' U1 h5 e

write.table(vgpm, vgpm.201603.txt, sep = \t, row.names = FALSE, quote = FALSE)

V' h$ w( W; A( k7 P

" L# M- k9 j& l6 @4 r s

转换后的数据框包括三列,分别是经度(lon)、纬度(lat)以及当前日期内该经纬度海区的NPPvar,单位mg C m-2 d-1)。

. D: c( K7 H$ e: i3 X7 E" ?0 x

使用nppr包匹配目标经纬度的遥感数据

& X# c( X3 [' b

默认情况下,下载的遥感数据是全球海洋的。nppr包同样提供了相关函数,便于我们从中提取特定区域的遥感数据,如下所示。

9 \ }3 k/ E* f; L

#获取指定日期和经纬度的遥感,例如在上述2016年3月的月平均NPP数据中提取120°E、20°N的NPP

% o% H( H3 q+ j z' @

match_sig(file.path = yourfolder, lon = 120, lat = 20, date = 2016-03-01)

6 ]; H7 m2 L9 r; w% ?: J# U W( f

#或者同时指定多个数据,不再多说

3 w/ f. U) y# u) b& u& Q4 M$ D- Z; C3 w

mydat <- data.frame(

: i8 ^6 g- s5 R( u9 K4 N+ H

lon = c(120, 112, 116),

8 C( \; p" Y p: F( m2 Y5 R4 h

lat = c(17, 15, 18),

( U O' m! r" c7 N4 d9 F4 ]+ B

date = c(2016-03-04, 2016-03-07, 2016-02-04)

1 U# E, Q! _6 t" L: U% i

)

B4 |0 a0 C0 k M# S3 L

match_df(mydat, file.path = yourfolder)

( d- [" t/ o+ H2 r8 d% v

绘制遥感地图

: e+ C! y( B( C

nppr包的函数geom_oce()可以用来绘制地图,例如我们作图展示来自遥感反演的NPP分布。

0 Z( @! k! {& D

#上述已经将下载的2016年3月的月平均NPP数据转换为数据框格式

; n8 o) s1 s% l* R

#我们仍以该数据为例作图,展示中国南海2016年3月的月平均NPP

/ ~# w% o$ U, }! F

library(viridis)

& Z3 U, e( P7 d: o6 S R' e+ M" Y2 I

library(ggplot2)

7 ~7 D% e' V b) E) v: I, s, j

ggplot()+

1 H4 M; h8 b+ ~" c

geom_oce(vgpm, aes(x = lon, y = lat, fill = var), lonlim = c(100, 120), latlim = c(7, 25))+

, _& f9 X! _' Q: `( i( y3 `# y- j8 e

scale_fill_viridis(option = D, direction = -1,breaks = seq(50, 1050, 100), limits = c(50, 1050))+

4 |: v1 G, Y' h& o9 u# b5 Q b

labs(x = Longitude, y = Latitude, fill = expression(NPP*~(*mg~C~m^-2~d^-1*)))

. }2 ^$ B! y/ A Z

- g5 C7 f. R- u8 r& }

根据时间和经纬度列表匹配遥感数据的批处理

- Y; _" c( F% X. \0 e" X4 m- }

实际情况中,经常需要对来自不同时间不同经纬度的大量站位匹配遥感,以下提供了一个批处理(不过这是自己先前瞎写的,然后一直偷懒一直用,俺也不知道写的对不对......写在这里仅为方便自己复制粘贴,大家慎用......

. r3 e( x* j0 a1 @! u$ q3 M

将待匹配的站位的经纬度、日期信息整合在一个文档中,如下所示的这样(本示例命名为“data.txt”)。

) y% F) F3 G6 I$ x$ |/ F. t

8 S3 X- V. ^8 }& M2 \

随后在R中读取该文件,设置一个循环,依次读取日期信息以下载当前日期的遥感(如月平均或8天平均的SSTPARChlaNPP等)。并再根据各站位的经纬度,从中匹配该站位附近的数据(比方说以0.1°为网格进行匹配,并将网格内的数据平均)。

6 j& J$ y# r+ Y5 c) K1 f

##如下以匹配SST数据为例做个演示

7 ~! M% r% W2 B

dat <- read.delim(data.txt)

8 I- L( D( f4 t: j

Date <- unique(dat$Date) #获取日期

; ]# R, @! P8 O

yourfolder <- paste0(getwd(), /, SST) #在当前工作路径下创建新路径以存放遥感数据

2 i/ q3 z* u8 q

dir.create(yourfolder)

& _/ j* ], [- L: H/ t

#通过循环依次获取各日期下的遥感(本示例以下载8天平均SST为例)

2 D( i: O! x) t. K' n, d

for (i in Date) {

9 v8 {/ N! ]4 B6 s1 B& y( H' B

yourfolder <- paste(getwd(), SST, i, sep = /)

3 \ T: S4 @9 w1 a: C: k: d

dir.create(yourfolder)

0 ]! B+ e1 A+ X! J5 ^( S6 P

get_sst(file.path = yourfolder, grid.size = low, time.span = dayly, satellite = MODIS, mindate = i, maxdate = i)

% O; h) g2 W3 I- r R$ q# ]7 R# w

yourfile <- dir(yourfolder)

: M: f( E# F% ? x- C4 g, b# H

hdf <- read_hdf(file.path = paste(yourfolder, yourfile, sep = /))

1 g& L5 n* W6 ^8 `1 r

write.table(hdf, paste(yourfolder, /, i, .xls, sep = ), sep = \t, row.names = FALSE, quote = FALSE)

! t( p0 A' H& R: \0 r8 i

}

" a& E B; W! I8 ?# b

#再根据列表中各站位的经纬度匹配当前日期的遥感(本示例计算0.1°网格内的平均)

: v) D* W. P) V" X. ]5 Z, @

for (i in 1:nrow(dat)) {

% T5 Y k) Z' |

Date <- dat[i,Date]

# Z$ e) r, m6 a$ }

yourfile <- paste(getwd(), /, SST, /, Date, /, Date, .xls, sep = )

; T4 c4 Q f5 w. A

hdf <- read.delim(yourfile)

. |8 D; c* z; D2 T( r

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), ]

9 [( `5 \9 V( h

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), ]

8 z: ~: ~' \6 ~8 x; t. g( {

dat[i,SST] <- mean(hdf$var)

/ y# \5 U+ o' _1 ?; B/ H, _. K

}

2 G! N4 y# ^+ N1 Y J4 c3 x# A

write.table(dat, SST+0.1.xls, sep =\t, quote = FALSE, row.names = FALSE)

, x. C( O$ h2 Q

2 g' l% |$ ]$ }) p9 Y& I# D

输出列表的最后一列添加了匹配的遥感数据(本示例匹配了SST)。

& T# U& V4 f! P1 L2 @/ P% }- S 8 Z. e1 f* q( I5 r+ O" R; r. G# V6 k $ Y: E' Y3 i" u. { % }# e0 e# I- A8 @* ^( \3 e
回复

举报 使用道具

全部回帖
暂无回帖,快来参与回复吧
懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
黑泽逢世
活跃在2025-1-26
快速回复 返回顶部 返回列表