CODEWITHSUNDEEP
rexcelextractnetcdfr-raster
marie_r
marie_r

Reputation: 67

Extract time series from netCDF in R

I created this file using TRMM_3B42_Daily product over 1998-01-01 to 1998-12-31. This is the script I used in R:

lon=seq(-91.875,-86.875,by= 0.25)
lat=seq(13.875,16.875,by= 0.25)

x_dim <- ncdim_def( "lon", "degrees_east", lon, create_dimvar=TRUE)
y_dim <- ncdim_def( "lat", "degrees_north", lat, create_dimvar=TRUE)
t_dim <- ncdim_def( "time", "days since 1997-12-31 12:00:00.0 -0:00", 1:365, unlim=FALSE)
mv=9999.900390625 
precipitation_var <- ncvar_def("precipitation", "mm", list(y_dim,x_dim,t_dim), mv)


nrow = 13 
ncol = 21 

NA.matrix=matrix(rep(NA,nrow*ncol)) 

precip=array(NA.matrix,c(nrow,ncol, 1))
for (i in 1:length(test01)){precip_nc=nc_open(test01[i])
precip_get_nc=ncvar_get(precip_nc,"precipitation") 
precip=abind(precip,precip_get_nc)}

precip=precip[,,-1]  

PRECIPITATION_nc = nc_create("PRECIPITATION_1998.nc", precipitation_var)

precipitation_nc_put=ncvar_put (PRECIPITATION_nc, precipitation_var, precip)

nc_close(PRECIPITATION_nc)

Following this link I tried extracting the values in order to plot a time series but it seems I am averaging the values of two cells instead of just extracting the values of a single cell. How do I fix this? Is there a way to create a loop so that it extracts the values of different cells? (in this case it would be 13 x 21 = 273)

b <- brick('PRECIPITATION_1998.nc')
be <- crop(b, extent(13.875, 14.125, -91.875,-91.625))
a <- aggregate(be, dim(be)[2:1], na.rm=TRUE)
v <- values(a)
write.csv(v, 'precip.csv', row.names=FALSE)

Also two other problems I found where that the dates in the excel file have an X in front and that the values are shown horizontally instead of vertically. Any help would be greatly appreciated!! Thanks

Upvotes: 2

Views: 2977

Answers (1)

lbusett
lbusett

Reputation: 5932

extraction of points data can be easily accomplished by creating a SpatialPoints object containing the point from which you want to extract data, followed by an extract operation. Concerining the other topics: The "X"s are added because column names can not start with numerals, so a character is added. The horizontal ordering can be easily changed after extraction with some transposing

This, for example, should work (It solves also the "X"s problem and changes the format to "column like"): 

library(raster)
library(stringr)
library(lubridate)
library(tidyverse)

b <- brick('/home/lb/Temp/buttami/PRECIPITATION_1998.nc')
lon = c(-91.875,-91.625)  # Array of x coordinates
lat <- c(13.875, 14.125)  # Array of y coordinates
points <- SpatialPoints(cbind(lat,lon)), # Build a spPoints object

# Etract and tidy
points_data <- b %>% 
  raster::extract(points, df = T) %>% 
  gather(date, value, -ID) %>% 
  spread(ID, value) %>%   # Can be skipped if you want a "long" table
  mutate(date = ymd(str_sub(names(b),2))) %>% 
  as_tibble()

points_data 

# A tibble: 365 × 3
         date   `1`   `2`
       <date> <dbl> <dbl>
1  1998-01-01     0     0
2  1998-01-02     0     0
3  1998-01-03     0     0
4  1998-01-04     0     0
5  1998-01-05     0     0
6  1998-01-06     0     0
7  1998-01-07     0     0
8  1998-01-08     0     0
9  1998-01-09     0     0
10 1998-01-10     0     0
# ... with 355 more rows

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

Upvotes: 2

Related Questions