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

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

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

9 K0 S. K& R, }0 ]

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

安装nppr包

# e$ K7 T' H, `: ]

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

5 `+ C4 y; v: P0 B8 L6 ~" a

#下载nppr包

: k- i6 ^. X" W* x) V/ k N

#install.packages(remotes)

! k& c# Z* p5 S6 v1 W& Y

remotes::install_github(chaoxv/nppr)

#加载相关R包

library(nppr)

library(RCurl)

library(XML)

library(R.utils)

8 A8 X7 }# J/ X8 a

library(tidyverse)

8 d" h& U# m( k

library(lubridate)

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

! {8 K/ n+ H8 S% j

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

; T" m' A' c- i# L, G2 O; H6 a5 t

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

) `8 T# W: K8 T- S G- w

#创建工作路径

yourfolder <- F:/R/nppr/vgpm

8 p) U6 v. C1 p5 D, L0 n

dir.create(yourfolder)

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

get_npp_vgpm(

file.path = yourfolder,

+ h! r# s4 P/ a3 @' y, f3 q

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

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

% s0 H$ \* k ]; e0 p) p" Z1 z* _& [2 U( u

satellite = MODIS, #选择卫星

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

maxdate = 2016-03-15

5 O! ^+ O! K7 Y3 g1 A4 T

)

在这个示例中，我们使用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数据做个转换

; X* D% W7 J9 ]2 o3 V$ f. j

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

vgpm <- read_hdf(file.path = yourfile)

head(vgpm)

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

$ X) y2 J$ O: P6 m$ |/ a- R

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

$ {; L2 M$ e% L, s3 g7 l6 s7 f) l: ?' Q

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

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

* R6 K4 H' x7 J2 {! t

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

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

* L; k) ]& p* N9 `

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

0 v3 h9 v% E5 O+ s

mydat <- data.frame(

lon = c(120, 112, 116),

4 l7 s( d8 I8 {9 t# M

lat = c(17, 15, 18),

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

)

. t/ F' L8 x1 D( V' R

match_df(mydat, file.path = yourfolder)

) f2 P9 m+ q; C6 i3 J1 m# D

绘制遥感地图

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

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

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

library(viridis)

9 @. B1 J# h- _: G4 Z- N! N6 ]

library(ggplot2)

! n$ O* Q/ Z- i4 h0 N( E& P* f7 a

ggplot()+

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))+

/ m9 G7 t* h( t' T. K& e$ P$ u" N' a

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

7 V) x1 X$ a/ t5 R

6 @# m' K9 t9 s- z4 k

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

: T- Q' {" t9 {: _# I4 }/ g/ t# v

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

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

: W! ~( p' f9 f4 J; S

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

# J/ @; K( Y. s

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

( l- R& n3 J2 {8 d0 m

dat <- read.delim(data.txt)

" m! L' {7 W. |) Y

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

m* q/ Y/ g4 @: H% U

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

+ c+ I% Y5 G/ x' r0 G7 E

dir.create(yourfolder)

+ A5 n+ _; l, f9 S4 U( i& s; C& P6 e1 _5 b

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

for (i in Date) {

& d9 Y/ k. G6 m8 @$ Z I! {

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

dir.create(yourfolder)

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

yourfile <- dir(yourfolder)

7 u8 I0 Q, t2 x2 X2 V

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

' j E5 X# e; n

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

O& p% }; u, x0 F! h

}

/ d H% p. |$ o) k0 b2 v, ^

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

& J, e) _* \* r q* u, y1 V2 p7 t

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

Date <- dat[i,Date]

( K" L3 Q+ m) |# |! k

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

! \2 Z; R$ @( v2 b5 K

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

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

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

}

6 p* }9 j& E$ J) K

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

/ I% a d4 V& J. Z% x0 ]

" q8 Y- V, D6 G" _4 c" @4 n a; k

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

" ^7 E1 R+ h+ q 5 K: a* O) p2 _/ a