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

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

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

2 C' @/ o$ T: `$ c5 F$ K. V

* S& |3 N8 Q& p/ x

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

, f6 d( Y; g8 r e4 `9 } t8 |. M

安装nppr包

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

* Y! o+ y- g0 ^, ?1 o6 ?- F2 f

#下载nppr包

#install.packages(remotes)

remotes::install_github(chaoxv/nppr)

* [; o( d* m5 P6 s

#加载相关R包

% k; X2 a/ ?4 |4 h

library(nppr)

library(RCurl)

library(XML)

5 \! j# V: |5 V6 b7 N

library(R.utils)

2 V$ e; w: K+ B/ ^ P+ i$ p

library(tidyverse)

+ S! x$ x! E. P6 ?4 u, o

library(lubridate)

& G: b7 @, o1 Y: Y( x

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

5 `- Y+ D& H. o3 |- @ n

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

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

5 t5 C5 U, B3 p0 y, M2 o2 h2 l) B

#创建工作路径

yourfolder <- F:/R/nppr/vgpm

dir.create(yourfolder)

' K$ f) m4 M6 w/ n

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

2 p# c7 F' g( w, K! A$ O1 u

get_npp_vgpm(

) C2 o" M8 J8 K; i4 R- E) V, C

file.path = yourfolder,

6 L! Z! H7 @! @: j$ V, s! z" r

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

2 i D* h0 p: I1 U& C3 @

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

/ Y* L" }3 D# N/ `7 H

satellite = MODIS, #选择卫星

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

% _# F# ~$ `/ @

maxdate = 2016-03-15

)

) e9 K! Q6 z2 t6 D) b

& E2 F3 q3 q% 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数据。

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

& u# _! B t9 v6 v

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

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

, a. j3 r5 w# U

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

vgpm <- read_hdf(file.path = yourfile)

head(vgpm)

$ e- k8 h- ^- V

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

$ k& O' F: b2 ]( L3 p1 g2 _

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

S% c) h5 a7 n8 \0 p7 R# m

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

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

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

6 i9 n6 Z5 Y7 A3 Z: X$ N8 s+ l

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

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

mydat <- data.frame(

lon = c(120, 112, 116),

7 w% s( b" M: B) k: V

lat = c(17, 15, 18),

% {1 P* j( S6 ]& z/ e% T" ?! g

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

0 _! n) {0 M* ]5 E9 t j( V% R

)

$ E" T' h+ j# x7 P$ f8 A

match_df(mydat, file.path = yourfolder)

绘制遥感地图

/ l' \2 I- S0 Z' c6 C% ?4 B+ e

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

, Y" g6 d$ u: s$ Z$ O

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

6 L$ a5 P- \3 w& ]* k

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

library(viridis)

" E# \1 h X' z* P- g+ Y

library(ggplot2)

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

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

* h! s l& h* O

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

' o, u7 W3 X7 k3 K7 t# v

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

7 U3 w$ x' k& g C: A4 P6 b, J

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

) E9 Y3 |- ~7 A1 f- @

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

' E0 Y! g! H& L% a2 |6 f3 K' f p

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

dat <- read.delim(data.txt)

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

* \ ]9 u- {3 L) u4 e6 ?6 A* b- Q

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

. F7 b) p: ]4 k+ c: I( h

dir.create(yourfolder)

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

3 I% O. L# y# \1 x5 z/ p# G& l, D

for (i in Date) {

# B0 ?: j+ s; p& M! `" D

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)

2 P1 v n# F- B* k0 q

yourfile <- dir(yourfolder)

4 g& |9 w0 t7 O( m; q+ L5 ~

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

7 {# N7 m* k- o# W0 m8 Y

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

}

5 `) o! t( C* L0 M

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

8 K; S5 H3 v5 a$ ?/ \) a- N6 j

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

Date <- dat[i,Date]

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

hdf <- read.delim(yourfile)

; N2 l' I; E- b/ H8 D/ 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), ]

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)

& W0 g% y% W( s& V7 j2 \

}

$ W( Y5 a' J# G

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

# m0 M6 z8 [: h

, a; F0 P3 {& t4 F4 L+ S' S

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

4 V9 g; K7 e2 r ; b% K# r. u) V