收藏本站 劰载中...网站公告 | 吾爱海洋论坛交流QQ群:835383472

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

[复制链接]
使用nppr包下载和处理海洋遥感数据(SST、Chla、PAR、NPP)- ^& s$ \! k, G$ d7 X1 u

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

f' M% _8 b6 c4 J' {. j% z3 p, i1 t5 @

, S" w$ U$ {3 D# u' G7 k f

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

( o3 H9 e' Q6 K! }+ R7 x R8 k

安装nppr包

6 N. L% H6 W* k

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

5 _7 g9 a4 \# `; E/ `7 ~

#下载nppr包

( {2 X) C3 y9 K5 [

#install.packages(remotes)

' E O, Q8 W7 ]" j' R

remotes::install_github(chaoxv/nppr)

) ?) a* \3 m7 N# I7 h2 W

#加载相关R包

0 O& B' p1 N4 c2 U

library(nppr)

4 e! K! ~; x" J' T9 e/ S

library(RCurl)

9 q5 T# H! X- A& G* i

library(XML)

3 c. S& W6 x0 X0 `

library(R.utils)

9 }9 `8 {8 W2 o; Z7 J

library(tidyverse)

: h/ V: D! C8 }( S; q

library(lubridate)

2 d4 \5 c, r! z

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

( a1 K5 U# ~# G. A# }! [

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

/ Y7 ?5 K x2 m

0 J- s& b) m% e4 h3 h% M

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

F+ n- p( X. l& F$ T( d+ O6 I

#创建工作路径

6 ]+ q1 G! z1 d8 \1 [6 Y x9 |

yourfolder <- F:/R/nppr/vgpm

$ F% G. u6 |. F9 Q

dir.create(yourfolder)

( o2 U6 Y3 @& h. u+ Y+ x- l

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

1 O# G6 B: c/ H I f5 S$ r G7 |

get_npp_vgpm(

& ^7 G6 e/ w# E* }+ H" s

file.path = yourfolder,

. o% l" A2 Q% q

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

2 Z1 ] z$ A; y6 g' M* X

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

! F X4 K; T0 {! M9 ~5 f% d5 a

satellite = MODIS, #选择卫星

) U( c. m. ^" s0 ?4 ?5 ?( }

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

@8 _0 P! N8 w8 ]5 N9 s/ F

maxdate = 2016-03-15

6 y: ^7 X# E! A% d. ~0 _

)

+ D7 L# t; j! |, r6 ]

3 L! d3 p2 M$ E @# ] h( {

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

- |, Y' D: C- s7 l1 `

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

3 n4 k2 @% {1 v- q: w7 J

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

9 K3 `9 a7 j' `0 y6 {1 p c( {4 T

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

8 A: x0 X0 A! }; D8 o

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

6 O. ^' u2 Z+ t9 S( g# H

vgpm <- read_hdf(file.path = yourfile)

& M5 P2 u B! b. j

head(vgpm)

! O. q8 f0 ?6 L' B* e, v

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

* @! G& u | w; P- f0 K5 c9 _+ @

# N ~4 i( n8 u6 I& e! W% [

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

: b1 {( s9 x- k e% o3 ]: l6 l

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

- J8 d+ Q6 g8 |# v9 T3 p

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

! g2 z' B3 ` j* p0 V; w. f

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

$ u: r; S; s [( l

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

, p# x( z9 r/ Q4 j9 B/ _

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

5 F* r% i* q, w

mydat <- data.frame(

8 J% W9 r' z' }

lon = c(120, 112, 116),

0 h4 I# o+ A$ q& S4 A8 E

lat = c(17, 15, 18),

" W# w* g! [. N" u; [7 u1 q

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

6 O: n; `) P0 [8 I

)

7 ^" Q4 f7 G: ] ~1 J& l

match_df(mydat, file.path = yourfolder)

+ g; w" A; k) j1 \' K# z+ m

绘制遥感地图

" x, q6 B$ K& R: T' e

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

" s1 H5 ?1 u; t9 L

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

$ t( R* |! x& ?8 R2 \- b. c

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

' a, Z, r+ H, E/ v# F; |

library(viridis)

' R1 B5 w. N+ ]" N) d+ B

library(ggplot2)

/ K) l: \+ S6 A; b9 x4 w

ggplot()+

( Q- b8 T4 ~0 [& h

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

9 N8 s' |& J) {; N# ~/ @

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

& v6 p; C$ _8 R% E' @ @

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

! j8 g' W* ^& S- |9 N

" c1 }4 f8 b( q+ E% A L- E7 @

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

; u2 y' T0 H: p2 K, j! N+ t, o

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

( ~: x2 g% f; ~; M: b. L. Y0 O( T

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

/ m9 b5 q* D; C$ ?

. E; n; m3 F9 Q' U5 J/ g/ n

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

+ U% _6 [* }# D6 l

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

' i0 l7 _2 d8 H7 M5 a& ]5 c

dat <- read.delim(data.txt)

% v) t& O: \/ F6 ]+ [* C( N; r' y

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

! y6 H/ [# k3 `' M7 N* w

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

, h4 e/ v& I7 Q0 S3 G! B8 p/ H

dir.create(yourfolder)

& Y0 q$ r+ Y8 d1 O1 K( n

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

- C- \% L" n5 P D

for (i in Date) {

" `" m# e: ?; f

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

* q( ], }* P+ \: M9 ?

dir.create(yourfolder)

. O4 g4 X) k( V$ w

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

% |& j9 T, k% X6 ?5 `7 |

yourfile <- dir(yourfolder)

6 f$ x+ r) k1 h" p* d9 o7 q

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

6 n' R( v: L$ Z: I5 G, k' z

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

( ?" h& T( I+ x8 @

}

- H) x( X* z: f8 ^

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

: \1 J" N7 X+ v0 i$ \1 D

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

1 h! K+ H$ T+ T

Date <- dat[i,Date]

# x/ m0 c+ T, @9 f# y6 }; g- j

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

$ y9 M( Y: j3 q7 A

hdf <- read.delim(yourfile)

' r& B4 N, ?4 X t T' 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 R: e* D3 t& g

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

) O0 k& K. {$ d4 d

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

8 f; u V+ N" B

}

& A( e* T2 G2 g3 I2 r; A

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

4 N1 F! t+ @& w& W- g9 V

" a% {5 V+ C' m- F% l6 ^1 b' I

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

* b4 N: D: n+ u. D4 ]) k$ V" E1 W' L& H9 h# t, G ! ^# n1 x9 T+ m- r( `6 C `. V$ k! o8 R; E$ p5 M/ o4 n: P7 v3 }( Y! C
回复

举报 使用道具

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