Scroll through the netcdf files and do the calculations - Python or R

Question

This is my first time using netCDF and I'm trying to wrap my head around working with it.

I have several netcdf version 3 files (NOAA NARR air.2m daily average values ​​all year long). Each file covers a year between 1979 and 2012. This is a 349 x 277 grid with a resolution of 32 km. Data has been downloaded from here .

Measurement is time (hours from 1/1/1800), and my variable of interest is air. I need to calculate the accumulated days with a temperature of 0. For example

Day 1 = +4 degrees, accumulated days = 0 Day 2 = -1 degrees, accumulated days = 1 Day 3 = -2 degrees, accumulated days = 2 Day 4 = -4 degrees, accumulated days = 3 Day 5 = +2 degrees, accumulated days = 0 Day 6 = -3 degrees, accumulated days = 1 

I need to save this data in a new netcdf file. I am familiar with Python and somewhat with R. What is the best way to cycle every day, check the value of previous days and based on this, output the value to a new netcdf file with exactly the same dimension and variable ... or just add another variable to the original netcdf file with the result I'm looking for.

Is it better to leave all files separate or merge them? I combined them with ncrcat and it worked fine, but the file is 2.3gb.

Thanks for the input.

My current progress in python:

 import numpy import netCDF4 #Change my working DIR f = netCDF4.Dataset('air7912.nc', 'r') for a in f.variables: print(a) #output = lat long x y Lambert_Conformal time time_bnds air f.variables['air'][1, 1, 1] #Output 298.37473 

To help me better understand what data structure I'm working with? Is ['air'] key in the above example and [1,1,1] are also keys? to get the value 298.37473. How can I then skip [1,1,1]?

Answer 1

You can use the very nice MFDataset function in netCDF4 to process multiple files as a single aggregated file without using ncrcat . Thus, the code will look like this:

 from pylab import * import netCDF4 f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc') # print variables f.variables.keys() atemp = f.variables['air'] print atemp ntimes, ny, nx = shape(atemp) cold_days = zeros((ny,nx),dtype=int) for i in xrange(ntimes): cold_days += atemp[i,:,:].data-273.15 < 0 pcolormesh(cold_days) colorbar() 

And here is one way to write a file (there may be simpler ways):

 # create NetCDF file nco = netCDF4.Dataset('/usgs/data2/notebook/cold_days.nc','w',clobber=True) nco.createDimension('x',nx) nco.createDimension('y',ny) cold_days_v = nco.createVariable('cold_days', 'i4', ( 'y', 'x')) cold_days_v.units='days' cold_days_v.long_name='total number of days below 0 degC' cold_days_v.grid_mapping = 'Lambert_Conformal' lono = nco.createVariable('lon','f4',('y','x')) lato = nco.createVariable('lat','f4',('y','x')) xo = nco.createVariable('x','f4',('x')) yo = nco.createVariable('y','f4',('y')) lco = nco.createVariable('Lambert_Conformal','i4') # copy all the variable attributes from original file for var in ['lon','lat','x','y','Lambert_Conformal']: for att in f.variables[var].ncattrs(): setattr(nco.variables[var],att,getattr(f.variables[var],att)) # copy variable data for lon,lat,x and y lono[:]=f.variables['lon'][:] lato[:]=f.variables['lat'][:] xo[:]=f.variables['x'][:] yo[:]=f.variables['y'][:] # write the cold_days data cold_days_v[:,:]=cold_days # copy Global attributes from original file for att in f.ncattrs(): setattr(nco,att,getattr(f,att)) nco.Conventions='CF-1.6' nco.close() 

If I try to view the resulting file in the Unidata NetCDF-Java Tools-UI GUI , this looks fine: Also note that here I just uploaded two datasets for testing, so I used

 f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc') 

as an example. For all data you can use

 f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.????.nc') 

or

 f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.*.nc') 

Answer 2

Here is the solution of R

 infiles <- list.files("data", pattern = "nc", full.names = TRUE, include.dirs = TRUE) outfile <- "data/air.colddays.nc" library(raster) r <- raster::stack(infiles) r <- sum((r - 273.15) < 0) plot(r) 

Answer 3

I know that it is rather late for this stream since 2013, but I just want to indicate that the decision taken does not provide a solution to the question posed. It seems that the question requires the length of each continuous temperature period to be below zero (note the question that the counter resets if the temperature exceeds zero), which may be important for climatic applications (for example, for agriculture), while the decision made gives everything only the total number of days in a year when the temperature is below zero. If this is really what mkmitchell wants (it was accepted as an answer), then this can be done from the command line in cdo without worrying about NETCDF input / output:

  cdo timsum -lec,273.15 in.nc out.nc 

so the loopback script would be:

 files=`ls *.nc` # pick up all the netcdf files in a directory for file in $files ; do # I use 273.15 as from the question seems T is in Kelvin cdo timsum -lec,273.15 $file ${file%???}_numdays.nc done 

If you want the total over the entire time period you can, then cat instead of _numdays files are much smaller:

 cdo cat *_numdays.nc total.nc cdo timsum total.nc total_below_zero.nc 

But again, the question seems to require that the accumulated days for the event be different, but not be provided with the accepted answer.