Skip to content

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

4 messages · Roy Mendelssohn - NOAA Federal, Eeusha Nafi, Michael Sumner

Original
Eeusha Nafi
# 
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!!*
Roy Mendelssohn - NOAA Federal
# 
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
**********************
"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.
Eeusha Nafi
# 
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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Michael Sumner
    

    

      
      #
    

  


  

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

        

            

      
      
      
    

        
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.