Creating a time series from a netcdf file for specific coordinates

Question

I am currently trying to create a time series of precipitation data from the CRU global land precipitation data. The file I am working has 60 time slices, from 2011-2015.

My code so far is:

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

The following error appears when I try to run this:

'Error in plot.window(...) : need finite 'ylim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf'

As you can see from above I am trying to create a timeseries for a range of coordinates, if it is simpler to find an average of these points and then plot a timeseries I would but also having difficulties with this. This is my first time trying to netCDF files so I not overly confident with it, any suggestions of how to do this would be greatly appreciated!

I do have some coordinates with NA values for precipitation if that's what you mean? I have used another code here: http://geog.uoregon.edu/bartlein/courses/geog607/Rmd/netCDF_01.htm which gets rid of the NAs but I haven't been able to combine it with the above code. — lmwedge267, Aug 01 '17 at 19:09
But the error message is about x axis (date), not about precipitation. You can hold x observations without y value, but not y values without the dependent x value. But I'll look the data, I've been working with netcdf in r in the past. Do you have a download link? — aldo_tapia, Aug 01 '17 at 19:16
Oh right sorry I didn't understand that. Yes I do, here is the link: https://www.dropbox.com/s/7urxt2l25mvl62i/cru_ts4.00.2011.2015.pre.dat.nc?dl=0 . Thank you very much for having a look. — lmwedge267, Aug 01 '17 at 19:30

score 1 · Answer 1 · answered Aug 01 '17 at 19:56

The problem is really simple, the first selected coordinate never has values. Is in an area without valid pixels. Let me show you:

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)

Only one of both points is over a valid pixel:

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...

So, you can plot only one column:

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

Okay I understand now thank you! If I want to check other pixel values will I able to add coordinates values at this point of the code: 'lon<-c(26.25,27.75) lat<-c(-16.25,-15.25)' and then change the xlim and ylim extent? — lmwedge267, Aug 01 '17 at 20:07
You only have to change `lon<-c(26.25,27.75) lat<-c(-16.25,-15.25)`. The `xlim` and `ylim` are for plotting purpose only — aldo_tapia, Aug 01 '17 at 20:10

score 1 · Answer 2 · answered Aug 02 '17 at 12:34

I think it is even easier to do this from the command line in CDO and then read the resulting time series file into R.

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

This picks up the nearest gridcell to your specified coordinates X/Y (nn="nearest neighbour"). This way you don't have to worry about defining a "box" around your location and worrying about whether this box contains no cells or more than one cell etc.

You can now read in the resulting timeseries in R.

ps: If you do want to define an area and make the timeseries for the average over that area then you want to use "sellonlatbox" and "fldmean":

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

This ensure the changing cell size as a function of latitude is accounted for (which using a simple arithmetic mean in R will not do, and thus will be in error).

Creating a time series from a netcdf file for specific coordinates

2 Answers2