0

我目前正在尝试从 CRU 全球陆地降水数据创建降水数据的时间序列。我正在工作的文件有 60 个时间片，从 2011 年到 2015 年。

到目前为止，我的代码是：

b<-brick('/Volumes/LDRIVE2/Masters/Dissertation/Global Surface Water Explorer/UEA_rainfall/cru_ts4.00.2011.2015.pre.dat.nc/cru_ts4.00.2011.2015.pre.dat.nc')

lon<-c(26.25,27.75)
lat<-c(-16.25,-15.25)
points<-SpatialPoints(cbind(lat,lon))

points_data <- b

raster::extract(points, df = T) %>% 
  gather(date, value, -ID) %>% 
  spread(ID, value) %>%    
  mutate(date = ymd(str_sub(names(b),2))) %>% 
  as_tibble()

plot(points_data$date,points_data$`1`)

当我尝试运行它时出现以下错误：

'plot.window(...) 中的错误：需要有限的 'ylim' 值此外：警告消息：1：在 min(x) 中：min 没有非缺失参数；返回 Inf 2：在 max(x) 中：max 没有非缺失参数；返回-Inf'

正如您从上面看到的那样，我正在尝试为一系列坐标创建一个时间序列，如果找到这些点的平均值然后绘制一个时间序列更简单，我会这样做但也有困难。这是我第一次尝试 netCDF 文件，所以我对它并不过分自信，任何关于如何做到这一点的建议将不胜感激！

4

2 回答 2

1

问题很简单，第一个选择的坐标从来没有值。位于没有有效像素的区域中。我来给你展示：

library(raster)
library(dplyr)
library(ndcf4)
library(tidyr)
library(lubridate)
library(stringr)
library(tibble)

b<-brick('~/Downloads/cru_ts4.00.2011.2015.pre.dat.nc')

lon<-c(26.25,27.75)
lat<-c(-16.25,-15.25)
points<-SpatialPoints(cbind(lat,lon))

plot(b[[1]], xlim=extent(points)[c(1,2)],ylim=extent(points)[c(3,4)])
plot(points, add=T)

在此处输入图像描述

两个点中只有一个在有效像素上：

points_data <- b %>% raster::extract(points, df = T) %>% 
  gather(date, value, -ID) %>% 
  spread(ID, value) %>%    
  mutate(date = ymd(str_sub(names(b),2))) %>% 
  as_tibble()

points_data %>% glimpse()

## Observations: 60
## Variables: 3
## $ date <date> 2011-01-16, 2011-02-15, 2011-03-16, 2011-04-16, 2011-05-16, 2011-06-16, 2011-07-16, 2011...
## $ 1    <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N...
## $ 2    <dbl> 115.099998, 6.400000, 37.799999, 70.599998, 8.200000, 0.000000, 0.000000, 0.000000, 1.100...

因此，您只能绘制一列：

plot(points_data$date,points_data$`2`)

在此处输入图像描述

于 2017-08-01T19:56:42.457 回答
1

我认为从 CDO 中的命令行执行此操作然后将生成的时间序列文件读入 R 会更容易。

cdo remapnn,lon=X/lat=Y crufile.nc timeseries.nc

这会选择离您指定的坐标 X/Y (nn="nearest neighbour") 最近的网格单元。这样您就不必担心在您的位置周围定义一个“框”，也不必担心这个框是否包含任何单元格或多个单元格等。

您现在可以在 R 中读取生成的时间序列。

ps：如果您确实想定义一个区域并为该区域的平均值制作时间序列，那么您想使用“sellonlatbox”和“fldmean”：

cdo fldmean -sellonlatbox,lon1,lon2,lat1,lat2 crufile.nc areamean.nc

这确保了考虑到作为纬度函数的单元格大小的变化（在 R 中使用简单的算术平均值不会这样做，因此会出错）。

于 2017-08-02T12:34:36.233 回答