Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Extract time series of a point ( lon, lat) from netCDF in R

Tags:

r

extract

netcdf

I am relatively new on R. I am trying to get time series of different points ( lat, lon) of temperature data from a netCDF file. My sample data file is here and here is the small file . I have tried netCDF package and the code i have used so far

library(ncdf)
obsdata = open.ncdf("obs.nc")

print.ncdf(obsdata) 

obsdatadates = obsdata$dim$time$vals
obsdatadates = as.Date(obsdatadates,origin = '1950-01-01') 
obsdatadates
obsoutput = get.var.ncdf(obsdata, varid = 'tasmin', start = c(1,1,1),
                         count = c(1,1,22280))
dim(obsoutput)
datafinal=merge(obsdatadates,obsoutput)

Can anyone help me to get a dataframe of timeseries ( first column) and value of data in another for a particular points( lat, lon) of that data. In this case I am looking for time series ( 1950-01-01 to 2010-12-31 for which the data is ) for a particular lat lon point ( and repeat for many points of interests) and for given variable(in this case tasmin). Your help would be appreciated. Thank you, aseem

like image 813
Cirrus Avatar asked Dec 16 '13 21:12

Cirrus


3 Answers

Perhaps use the raster package, this won't work for all NetCDF files but it does for yours:

library(raster)
## brick reads all 22280 layers
r <- brick("obs.nc", varname = "tasmin")
## extract works for all time steps
vals <- extract(r, matrix(c(-120, 52.5), ncol = 2))

dim(vals)
## [1]     1 22280

Note that gives a 1-row, many column matrix because I only gave a single point to extract().

(The extraction is simple with direct copy from the nearest cell, use method = "bilinear" to do interpolation). See ?extract for other options.

like image 175
mdsumner Avatar answered Nov 14 '22 23:11

mdsumner


Here's how I would proceed to do this with ncdf:

library(ncdf)
obsdata = open.ncdf("obs1.nc")
obsdatadates = as.Date(obsdata$dim$time$vals,origin = '1950-01-01')
#Get the whole data first
obsoutput = get.var.ncdf(obsdata, varid = 'tasmin')
#Prepare your points of interest
points_of_interest = data.frame(lat=seq(1,8,by=2),lon=c(1,5,3,6))
#Subset your data accordingly
data_at_point = apply(points_of_interest,1,function(x)obsoutput[x[1],x[2],])
#Turn it into a dataframe
data_at_point = as.data.frame(data_at_point)
#Add the dates to the dataframe
data_at_point$Date = obsdatadates
like image 42
plannapus Avatar answered Nov 14 '22 23:11

plannapus


the 'ncdf'-package is deprecated: http://cirrus.ucsd.edu/~pierce/ncdf/

update with the ncdf4 package:

library(ncdf4)
obsdata <- nc_open("obs1.nc")
print(obsdata) # check that dims are lon-lat-time

# location of interest
lon <- 6  # longitude of location
lat <- 51 # latitude  of location

# get dates
obsdatadates <- as.Date(obsdata$dim$time$vals, origin = '1950-01-01')

# get values at location lonlat
obsoutput <- ncvar_get(obsdata, varid = 'tasmin',
                  start= c(which.min(abs(obsdata$dim$longitude$vals - lon)), # look for closest long
                           which.min(abs(obsdata$dim$latitude$vals - lat)),  # look for closest lat
                           1),
                  count = c(1,1,-1)) #count '-1' means 'all values along that dimension'that dimension'
# create dataframe
datafinal <- data.frame(dates= obsdatadates, obs = obsoutput)
like image 30
Lennert Avatar answered Nov 14 '22 21:11

Lennert