R help - Aug 2017 - Extracting subset from netCDF file using lat/lon and converting into .csv in R

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.

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]]

Two questions:

1.  Is the order of the dimensions shown what is shown if you look at str(ncin)
- I mean shown at the end where it describes the variable and its dimensions?

2.  Is you problem how to subset the netcdf file,  how to write to the .csv
file, or both?

-Roy
> On Aug 28, 2017, at 2:21 PM, 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.
> 
> 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.
> 
**********************
"The contents of this message do not reflect any position of the U.S.
Government or NOAA."
**********************
Roy Mendelssohn
Supervisory Operations Research Analyst
NOAA/NMFS
Environmental Research Division
Southwest Fisheries Science Center
***Note new street address***
110 McAllister Way
Santa Cruz, CA 95060
Phone: (831)-420-3666
Fax: (831) 420-3980
e-mail: Roy.Mendelssohn at noaa.gov www: http://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill."
"From those who have been given much, much will be expected" 
"the arc of the moral universe is long, but it bends toward justice"
-MLK Jr.

Thank you for the kind reply.
1. The order of the dimension is reflected when used as print(ncin)
function. This function has been used to have a quick look at the structure
of the default netCDF file.
2. Firstly, I am in need of subsetting this netCDF file and then have to
convert the slice into .csv format.

So far, I could read the netCDF file (variables and dimensions) and
fragment the time into fields. Therefore, I would be grateful to get
further help to reach the goal.

Regards,
Eeusha

On Mon, Aug 28, 2017 at 11:56 PM, Roy Mendelssohn - NOAA Federal <
roy.mendelssohn at noaa.gov> wrote:
> Two questions:
>
> 1.  Is the order of the dimensions shown what is shown if you look at
> str(ncin) - I mean shown at the end where it describes the variable and its
> dimensions?
>
> 2.  Is you problem how to subset the netcdf file,  how to write to the
> .csv file, or both?
>
> -Roy
>
> > On Aug 28, 2017, at 2:21 PM, 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.
> >
> > 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.
> >
>
> **********************
> "The contents of this message do not reflect any position of the U.S.
> Government or NOAA."
> **********************
> Roy Mendelssohn
> Supervisory Operations Research Analyst
> NOAA/NMFS
> Environmental Research Division
> Southwest Fisheries Science Center
> ***Note new street address***
> 110 McAllister Way
> Santa Cruz, CA 95060
> Phone: (831)-420-3666
> Fax: (831) 420-3980
> e-mail: Roy.Mendelssohn at noaa.gov www: http://www.pfeg.noaa.gov/
>
> "Old age and treachery will overcome youth and skill."
> "From those who have been given much, much will be expected"
> "the arc of the moral universe is long, but it bends toward
justice" -MLK
> Jr.
>
>
	[[alternative HTML version deleted]]

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]]