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

使用nppr包下载和处理海洋遥感数据（SST、Chla、PAR、NPP） S$ P3 E4 J, k {: A# _! m' n

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

$ i0 w: e/ P3 |6 Q# ]

$ c* t U7 i3 f- H6 ~) [% t5 u; t

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

$ |6 I$ w: c% B& \4 z: f% U+ Q

安装nppr包

可在github（https://github.com/chaoxv/nppr）获取nppr包。

" q6 ^) C1 s5 H8 g# O

#下载nppr包

#install.packages(remotes)

remotes::install_github(chaoxv/nppr)

& r. H, T; l3 b0 u* O. k9 H

#加载相关R包

2 M, L8 s8 \/ Y* H& G* w) K

library(nppr)

; n! S$ j8 ^. ]4 Y" T" C. n

library(RCurl)

library(XML)

Y/ b1 ]6 q% ~, y$ e; A% f4 j

library(R.utils)

library(tidyverse)

library(lubridate)

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

; G2 h Y/ J# b

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

$ D4 d p+ l- ?) T" L2 V, n

9 _ {& b% Z' [ V- h! |$ ]

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

9 ]- L* G1 G1 b9 \! b' [

#创建工作路径

8 w$ I6 y0 X/ `) i# \" `$ C

yourfolder <- F:/R/nppr/vgpm

dir.create(yourfolder)

#以VGPM（NPP的一种遥感算法）为例做个演示，详情?get_npp_vgpm

get_npp_vgpm(

- |' f6 z- j6 I7 T! q; } }

file.path = yourfolder,

6 Q" c' x4 B; ^4 ?# c& P( p# [$ v

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

! v. K0 `. D1 k; H* y( G+ h4 z

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

satellite = MODIS, #选择卫星

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

- n4 Y3 X$ ]1 i. E0 K

maxdate = 2016-03-15

" n; b% N# _, x

)

" i0 Y [9 b% t- }. M

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

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

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

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

K; i$ Z E4 v. u

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

' d i6 M5 {7 }' X2 d& N0 i

vgpm <- read_hdf(file.path = yourfile)

4 N2 M, J. T2 `$ ?2 v. J: t7 ^

head(vgpm)

" P' i* y/ I! j' J( I

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

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

1 Q7 f, p$ n' x

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

, p4 T: n- I* {* O

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

& g# h- G8 w0 v: X9 z

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

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

: P: S; H0 v1 k: D5 S1 v& q

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

mydat <- data.frame(

! j A9 h$ G. Z

lon = c(120, 112, 116),

7 Y9 c7 ^0 d1 M- I7 h

lat = c(17, 15, 18),

1 n. y9 t' T6 m" E2 S4 a4 g; {6 B" d

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

)

match_df(mydat, file.path = yourfolder)

绘制遥感地图

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

0 H( W% [0 x% W* o7 ~" z8 N

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

# o5 [* B7 s1 `9 l

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

library(viridis)

6 G8 x& {- L3 E( }% S/ b8 {

library(ggplot2)

ggplot()+

! [. [% D/ u8 s+ |5 D& O

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

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

- `+ a8 _* R% Q) U8 d3 y: P

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

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

- N( Z2 ^. D$ l9 ?. ~/ ~

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

9 j u. I8 ?: T+ k

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

. P* x. G+ F3 i. [

) Q3 q! d o, J5 h

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

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

* _- P8 D6 }) y/ y# \

dat <- read.delim(data.txt)

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

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

. ?$ R+ `: Y- @2 P) C

dir.create(yourfolder)

% r* o4 `3 ]' Y4 d) m' `

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

( Y+ d. q) z2 z5 u! q

for (i in Date) {

! L; g" @5 K# E4 S" X/ L- Z2 _

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

: n: ?0 Q' g) D# C0 d

dir.create(yourfolder)

6 S4 U- `* w% @- g3 b4 m

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

& i9 k+ f: ^& Y, p2 U. `: x$ }

yourfile <- dir(yourfolder)

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

7 a6 A& `4 }. }. l. E8 U& U; I& A

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

}

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

- k) q8 K8 G1 G; e A$ u4 T0 i

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

; P0 ~/ f" t1 L' @( p$ y! E

Date <- dat[i,Date]

, d; Y+ G d8 d8 |5 K5 C

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

hdf <- read.delim(yourfile)

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

% L+ U8 f, B e5 U$ {) l$ F3 z1 E: T

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

5 ^' v0 o6 w$ L1 m

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

}

) l1 I9 x8 p. t; F' s0 y

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

- g2 E& s1 o, h/ t' {5 h2 w. \

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

4 T$ l5 v$ `3 R4 p: N+ l0 ^. D $ Z, q* A' e$ V 4 ]: Z# T$ r/ F* | n' J: M Y8 Q