Extract time series of a point ( lon, lat) from netCDF in R
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.
Extract time series from netCDF in R
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`)
Use R to extract time series from netcdf data
This might do it, where "my.variable"
is the name of the variable you are interested in:
library(survival)
library(RNetCDF)
library(ncdf)
library(date)
setwd("c:/users/mmiller21/Netcdf")
my.data <- open.nc("my.netCDF.file.nc");
my.time <- var.get.nc(my.data, "time")
n.latitudes <- length(var.get.nc(my.data, "lat"))
n.longitudes <- length(var.get.nc(my.data, "lon"))
n.dates <- trunc(length(my.time))
n.dates
my.object <- var.get.nc(my.data, "my.variable")
my.array <- array(my.object, dim = c(n.latitudes, n.longitudes, n.dates))
my.array[,,1:5]
my.vector <- my.array[1, 2, 1:n.dates] # first latitude, second longitude
head(my.vector)
baseDate <- att.get.nc(my.data, "NC_GLOBAL", "base_date")
bdInt <- as.integer(baseDate[1])
year <- date.mdy(seq(mdy.date(1, 1, bdInt), by = 1,
length.out = length(my.vector)))$year
head(year)
Extracting values for specific lat long from netcdf
You can use the following code for data extraction from .nc files
dat <- structure(list(locatioID = paste0('ID', 1:16), lon = c(73.73, 86, 73.45, 86.41, 85.36, 81.95, 82.57, 75.66, 82.03,
81.73, 85.66, 85.31, 81.03, 81.70, 87.03, 73.38),
lat = c(24.59, 20.08, 22.61, 23.33, 23.99, 19.09, 18.85, 15.25, 26.78,
16.63, 25.98, 23.28, 24.5, 21.23, 25.08, 21.11)),
row.names = c(1L, 3L, 5L, 8L, 11L, 14L, 17L, 18L, 19L, 21L,
23L, 26L, 29L, 32L, 33L, 35L), class = "data.frame")
temp <- brick("chirps-v2.0.1981.days_p05.nc")
xy <- dat[,2:3] #Column 1 is longitude and column 2 is latitude
xy
spts <- SpatialPoints(xy, proj4string=CRS("+proj=longlat +datum=WGS84"))
#Extract data by spatial point
temp2 <- extract(temp, spts)
temp3 <- t(temp2) #transpose raster object
colnames(temp3) <- dat[,1] #It would be better if you have the location names corresponding to the points
head(temp3)
write.csv(temp3, "Rainfall.csv")
Extracting SST time series at multiple lat, lon locations using CDO
So here is my attempt to answer the question as it was stated in the comments (i.e. you wanted an index which was the midshelf locations averaged and then subtracting the same latitude SST sampled at Longitude=9E). I assume the locations are stored pair-wise in a text file called "locations.txt" as in your question above. The loop part of the answer is from one of this question's solutions.
# first loop over the pairs of indices in the text files.
while read -r -a fields; do
for ((i=0; i < ${#fields[@]}; i += 2)); do
# precise lon/lat for sampled mid-shelf
cdo remapnn,"lon=${fields[i]}/lat=${fields[i+1]}" in.nc pt_${i}.nc
# same lat but lon=9E (change as wanted)
cdo remapnn,"lon=9/lat=${fields[i+1]}" in.nc 9E_${i}.nc
done
done < example.txt
# now take the ensemble average over the points.
cdo ensmean pt_*.nc mid_shelf_sst.nc
cdo ensmean 9E_*.nc mid_shelf_9E.nc
# and calculate the index
cdo sub mid_shelf_sst.nc mid_shelf_9E.nc SST_index.nc
Related Topics
Ggplot2 Y Axis Label Decimal Precision
Programmatically Rename Columns in Dplyr
Multiple Filled.Contour Plots in One Graph Using with Par(Mfrow=C())
How Can Library() Accept Both Quoted and Unquoted Strings
Want Only the Time Portion of a Date-Time Object in R
Ggplot2 Error "No Layers in Plot"
R - How to Add Row Index to a Data Frame, Based on Combination of Factors
How to Plot Mean and Standard Error in Boxplot in R
R: Data.Table Count !Na Per Row
Fama MACbeth Standard Errors in R
Rolling Regression by Group in the Tidyverse
Mapping Specific States and Provinces in R
How to Jitter Two Ggplot Geoms in the Same Way
Oauth Authentification to Fitbit Using Httr
How to Find the First and Last Occurrences of an Element in a Data.Frame