[R] Extracting subset from netCDF file using lat/lon and converting into .csv in R

[Michael Sumner]

On Tue, 29 Aug 2017 at 07:22 Eeusha Nafi <eshad002 at gmail.com> wrote:

> I have a series of nertCDF files containing global data for a particular
> variable, e.g. tmin/tmax/precipiation/windspeed/relative
> humuidity/radiation etc. I get the following information when using
> *nc_open* function in R:
>
> datafile: https://www.dropbox.com/s/xpo7zklcmtm3g5r/gfdl_preci.nc?dl=0
>
> File gfdl_preci.nc (NC_FORMAT_NETCDF4_CLASSIC):
>
> 1 variables (excluding dimension variables):
>         float prAdjust[lon,lat,time]
>             _FillValue: 1.00000002004088e+20
>             missing_value: 1.00000002004088e+20
>             comment: includes all types (rain, snow, large-scale,
> convective, etc.)
>             long_name: Bias-Corrected Precipitation
>             units: kg m-2 s-1
>             standard_name: precipitation_flux
>
>      3 dimensions:
>         lon  Size:720
>             standard_name: longitude
>             long_name: longitude
>             units: degrees_east
>             axis: X
>         lat  Size:360
>             standard_name: latitude
>             long_name: latitude
>             units: degrees_north
>             axis: Y
>         time  Size:365   *** is unlimited ***
>             standard_name: time
>             units: days since 1860-1-1 00:00:00
>             calendar: standard
>             axis: T
>
>     14 global attributes:
>         CDI: Climate Data Interface version 1.7.0 (
> http://mpimet.mpg.de/cdi)
>         Conventions: CF-1.4
>         title: Model output climate of GFDL-ESM2M r1i1p1 Interpolated
> to 0.5 degree and bias corrected using observations from 1960 - 1999
> for EU WATCH project
>         CDO: Climate Data Operators version 1.7.0 (
> http://mpimet.mpg.de/cdo)
>         product_id: input
>         model_id: gfdl-esm2m
>         institute_id: PIK
>         experiment_id: historical
>         ensemble_id: r1i1p1
>         time_frequency: daily
>         creator: isimip at pik-potsdam.de
>         description: GFDL-ESM2M bias corrected impact model input
> prepared for ISIMIP2
>
> Now I want to extract a subset from this dataset using pair of lon and lat
> points, e.g., (12.875, -11.625) & (8.875, 4.125) and convert the file into
> .csv format.
>


For this file I would recommend using higher level tools, you need the
raster and ncdf4 packages. This data "maps" into the regular grid
assumptions required by raster, and it's pretty trivial. In the end you
really only need raster::brick and raster::extract to do this task, but
I've added a bit of context as it's probably all pretty unfamiliar.

f <- "~/gfdl_preci.nc"

library(raster)
pr <- brick(f)

## this is a multilayer raster in longlat
print(pr)

## our query points are in the same coordinate system
## of the brick, and they fall within the domain (assuming x,y, long,lat)
pt <- cbind(c(12.875, 8.875), c(-11.625, 4.125))
## check our assumptions with an actual map of the first slice
plot(pr[[1]])
points(pt)
## we have a nrow(pt) * nlayer(pr) array extracted
ex <- raster::extract(pr, pt)
dim(ex)
head(t(ex))
## note how one point falls out of the data range (in the sea?)
matplot(t(ex), type = "l")

Cheers, Mike.



>
> so far I could make up to this step:
>
> # load the ncdf4 package
> library(ncdf4)
> # set path and filename
> setwd("D:/netcdf")
> ncname <- "gfdl_preci"
> ncfname <- paste(ncname, ".nc", sep = "")
> dname <- "prAdjust"
> # open a netCDF file
> ncin <- nc_open(ncfname)
> print(ncin)# get longitude and latitude
> lon <- ncvar_get(ncin,"lon")
> nlon <- dim(lon)
> head(lon)
>
> lat <- ncvar_get(ncin,"lat")
> nlat <- dim(lat)
> head(lat)
>
> print(c(nlon,nlat))
> # get time
> time <- ncvar_get(ncin,"time")
> time
>
> tunits <- ncatt_get(ncin,"time","units")
> nt <- dim(time)
> nt
> tunits
> # get variable
> preci.array <- ncvar_get(ncin,dname)
>
> dlname <- ncatt_get(ncin,"prAdjust","long_name")
>
> dunits <- ncatt_get(ncin,"prAdjust","units")
>
> fillvalue <- ncatt_get(ncin,"prAdjust","_FillValue")
>
> dim(preci.array)
> # split the time units string into fields
> tustr <- strsplit(tunits$value, " ")
>
> tdstr <- strsplit(unlist(tustr)[3], "-")
>
> tmonth = as.integer(unlist(tdstr)[2])
>
> tday = as.integer(unlist(tdstr)[3])
>
> tyear = as.integer(unlist(tdstr)[1])
>
> chron(time, origin = c(tmonth, tday, tyear))
>
> *Any help would be appreciated!!*
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]