CODEWITHSUNDEEP
rtimenetcdfcdo-climate
lmwedge267
lmwedge267

Reputation: 1

Creating a time series from a netcdf file for specific coordinates

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!

Upvotes: 0

Views: 1652

Answers (2)

ClimateUnboxed
ClimateUnboxed

Reputation: 8057

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

Upvotes: 1

aldo_tapia
aldo_tapia

Reputation: 1228

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)

enter image description here

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

enter image description here

Upvotes: 1

Related Questions