Re: help creating a time netcdf file

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Mon Mar 17 2014 - 07:57:52 MDT

[1] The model is (likely) written in fortran.
     Different languages & formats use different ordering.
     Fortran is column major; netCDF (NCL also) is row major.

     fortran: x(mlon,nlat,ntim) <==> x(ntim,nlat,mlon) :netCDF [NCL]

     Different languages use different ordering.
     While they appear complete opposites, they refer to the
     same array memory layout. See [5]

[2] Did you do what I suggested? re: ncl_filedump -itime ....
     Please do so and *carefully look* at the output.
     Note the 'time' dimension.

[3] The unlimited dimension is needed by certain tools
     (eg: netCDF operators: ncrcat). NCL does not care and (likely) the
      model does not care.

[4] Please read:
http://www.ncl.ucar.edu/Document/Functions/Built-in/setfileoption.shtml
http://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml

[5]
********************************************************************
If you are new to NCL please read the Mini-Language Manual at:
        http://www.ncl.ucar.edu/Document/Manuals/

If you must reorder an array, please see
        Section 2.11.1 Dimension reordering:

Section 7.7 has an illustration of NCL/Fortran array mapping
*********************************************************************

[6] attached is a simple [untested] script to read multiple files,
     extract a sub-region, create netCDF with file attributes

On 3/15/14, 10:05 AM, john koudelka wrote:
> Dennis and Wei,
>
> Thanks for your responses.
>
> I'm working with a set of grib files, where each file is a step in a time.
> So, for example, my directory contains a file for each hour:
> NLDAS_FORA0125_H.A20101001.0000.002.grb
> NLDAS_FORA0125_H.A20101001.0100.002.grb
> NLDAS_FORA0125_H.A20101001.0200.002.grb
> . . .
> etc.
>
> Each NLDAS file contains a 2d temperature dataset (I know it contains other
> variables, temp is all I'm trying to get right now).
>
> I'm trying to loop through the directory, extract, clip, and regrid each
> temperature dataset, then place them into a single netcdf file. Ie,
> temp.ncwould contain each of the temperature grids, in (y,x,t) format.
>
> I don't need an unlimited dimension in time, since I obtain a file count at
> the beginning of the script, subsequently allowing me to declare the size
> of the time dimension.
>
> The model I'm working with is only compatible with netcdf 3, and prefers
> time located in the third dimension, hence the insistence to setup the file
> as y,x,t. From what I've read, there seems to be difficulty in setting it
> up like this instead of t, y, x. Is there a function that would permit me
> to swap dimensions around? Or is this only a limitation when trying to set
> the 3rd dimension as unlimited? Since I know the size of the time
> dimension, will I be OK?
>
> I know this is a bit involved, and I'm still learning NCL and netCDF
> concepts, so thank you for your patience and thoughtful guidance.
>
>
> John
>
>
>
>
>
> On Fri, Mar 14, 2014 at 4:25 PM, Dennis Shea <shea@ucar.edu> wrote:
>
>> In netCDF-3, the unlimited dimension must be the record
>> dimension. The record dimension is the leftmost dimension.
>> You could try commenting or deleting that or as Wei has suggested
>>
>> =======> Wei Huang
>> Just try to change line:
>>
>>
>> dimUnlim = (/False,False,True/)
>>
>> to
>>
>> dimUnlim = (/False,False,False/)
>>
>> Wei
>> ========
>>
>>
>> or
>>
>> It is my speculation that the GRIB file you are reading
>> has only one time step. By default, NCL will not read
>> degenerate (singleton) dimensions. Hence, GRIB files with
>> only one time step will not have an explicit 'time' dimension.
>>
>> For illustration ... Try
>>
>> http://www.ncl.ucar.edu/Document/Tools/ncl_filedump.shtml
>>
>> Specifically ....
>>
>> %> ncl_filedump sample_file.grb | less
>>
>> Then try with the -itime option
>>
>> %> ncl_filedump -itime sample_file.grb | less
>>
>> The latter forces NCL to create singleton dimensions.
>> ---
>>
>> Programmatically .... use setfileoption
>> http://www.ncl.ucar.edu/Document/Functions/Built-in/setfileoption.shtml
>>
>>
>> ; Force a time dimension for single element time
>> setfileoption("grb","SingleElementDimensions","Initial_time")
>>
>> ; open the ERA-I GRIB files from NCAR-CISL RDA
>> g = addfile("sample_file.grb" , "r")
>> t =g->TMP....
>> printVarSummary(t)
>>
>> or better
>>
>> latN = 42.852
>> latS = 40.500
>> lonL = -112.620
>> lonR = -110.480
>>
>> t =g->TMP....(:,{latS:latN},{lonL:lonR})
>> printVarSummary(t)
>>
>> you can rename dimensions as needed
>>
>> t!0 = "timE
>> t!1 = "lat"
>> t!2 = "lon"
>>
>> We would recommend writing the netCDF the 'simple way'
>>
>> http://www.ncl.ucar.edu/Applications/method_1.shtml
>>
>> ---
>> system("/bin/rm -f simple.nc") ; remove any pre-existing file
>> ncdf = addfile("simple.nc" ,"c") ; open output netCDF file
>>
>> ; make time and UNLIMITED dimension
>> ; recommended for most applications
>> filedimdef(ncdf,"time",-1,True)
>>
>> ; output variables directly
>> ncdf->T = t
>>
>>
>>
>>
>> On 3/14/14, 3:59 PM, john koudelka wrote:
>>
>>> Wei, thanks for the suggestions.
>>>
>>>
>>> The script errors with:
>>> fatal:NetCDF: NC_UNLIMITED in the wrong index
>>> fatal:["NclFile.c":474]:FileAddVar: an error occurred while adding a
>>> variable to a file, check to make sure data type is supported by the
>>> output
>>> format
>>>
>>> And refers to this line:
>>> filevardef(ncdf_out, "TEMP" ,typeof(tc) ,getvardims(tc))
>>>
>>>
>>>
>>>
>>>
>>> On Fri, Mar 14, 2014 at 3:41 PM, Wei Huang <huangwei@ucar.edu> wrote:
>>>
>>> John,
>>>>
>>>> Try the script attached below.
>>>>
>>>> Wei
>>>>
>>>> ---------
>>>>
>>>> ;**************************************************
>>>> begin
>>>> ncfiles=systemfunc("ls *.nc")
>>>> if(dimsizes(ncfiles) .gt. 0) then
>>>> system("rm *.nc") ;remove any pre-existing file
>>>> end if
>>>>
>>>> ;specify extent
>>>> latN=42.852
>>>> latS=40.500
>>>> lonL=-112.620
>>>> lonR=-110.480
>>>>
>>>> cmd="ls *.grb"
>>>> grbFiles=systemfunc(cmd)
>>>> nfiles=dimsizes(grbFiles)
>>>>
>>>> ;create time var
>>>> ;should be 'hours since 2010-10-01T00.00
>>>> ;time should be provided as floating numbers
>>>> start=str_split(grbFiles(0),".")
>>>> year=str_get_cols(start(1), 1, 4)
>>>> month=str_get_cols(start(1), 5, 6)
>>>> day=str_get_cols(start(1), 7, 8)
>>>> hour=str_get_cols(start(2),0,1)
>>>>
>>>> time = fspan(0,nfiles,1)
>>>> time!0 = "time"
>>>> time@units = "hours since "+year+"-"+month+"-"+day+"T00.00"
>>>> time&time = time
>>>> timec=tocharacter(time)
>>>> ;print(time)
>>>>
>>>>
>>>> ;create and define the output file
>>>> ncdf_out = addfile("./test.nc" ,"c")
>>>> ncdf_out@title="Temperature Data: degC"
>>>> ncdf_out@source="NLDAS Forcing2 - hourly"
>>>> ncdf_out@creation_date = systemfunc ("date")
>>>> setfileoption(ncdf_out,"DefineMode",True)
>>>>
>>>>
>>>> ;get dimensions and setup file
>>>> g=addfile(grbFiles(0),"r")
>>>> l=g->lat_110
>>>> lat=l({latS:latN})
>>>> lat!0="lat"
>>>> lon=g->lon_110({lonL:lonR})
>>>> lon!0="lon"
>>>> nlat=dimsizes(lat)
>>>> nlon=dimsizes(lon)
>>>>
>>>> dimNames = (/"lat","lon","time"/)
>>>> dimSizes = (/nlat,nlon,dimsizes(time)/)
>>>> dimUnlim = (/False,False,True/)
>>>> filedimdef(ncdf_out,dimNames,dimSizes,dimUnlim)
>>>>
>>>>
>>>> ;setup the variable
>>>> t=g->TMP_110_HTGL
>>>> tc = new((/nlat,nlon,dimsizes(time)/), float)
>>>> tc!0="lat"
>>>> tc!1="lon"
>>>> tc!2="time"
>>>> tc&time=time
>>>> tc&lat=lat
>>>> tc&lon=lon
>>>> tc@sub_center=t@sub_center
>>>> tc@center=t@center
>>>> tc@long_name=t@long_name
>>>> tc@units="degC"
>>>> tc@level_indicator=t@level_indicator
>>>> tc@gds_grid_type=t@gds_grid_type
>>>> tc@parameter_table_version=t@parameter_table_version
>>>> tc@parameter_number=t@parameter_number
>>>> tc@model=t@model
>>>> tc@level=t@level
>>>>
>>>> ;setup file variables
>>>> filevardef(ncdf_out, "lat", typeof(lat), getvardims(lat))
>>>>
>>>> filevardef(ncdf_out, "lon", typeof(lon), getvardims(lon))
>>>> filevardef(ncdf_out, "time", typeof(time),getvardims(time))
>>>>
>>>> filevardef(ncdf_out, "TEMP", typeof(tc), getvardims(tc))
>>>>
>>>> ncdf_out->lat = (/lat/)
>>>> ncdf_out->lon = (/lon/)
>>>> ncdf_out->time = (/time/)
>>>>
>>>> ; loop through the files and process
>>>> do n = 0, nfiles-1
>>>> g=addfile(grbFiles(n),"r")
>>>>
>>>> ;do the conversions and add to ncdf_out
>>>> ; TEMPERATURE
>>>> t=g->TMP_110_HTGL
>>>> tc=t({latS:latN},{lonL:lonR})-273.15
>>>> ;tc(,,n)=t({latS:latN},{lonL:lonR})-273.15
>>>> ;this spot reserved for regridding to higher resolution
>>>> ncdf_out->TEMP(:,:,n)=(/tc/)
>>>> end do
>>>> print(ncdf_out)
>>>> end
>>>>
>>>>
>>>> huangwei@ucar.edu
>>>> VETS/CISL
>>>> National Center for Atmospheric Research
>>>> P.O. Box 3000 (1850 Table Mesa Dr.)
>>>> Boulder, CO 80307-3000 USA
>>>> (303) 497-8924
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Mar 14, 2014, at 2:15 PM, john koudelka <john.koudelka@gmail.com>
>>>> wrote:
>>>>
>>>> I need some guidance as to how to properly consolidate a time series of
>>>> temperature grids into a single netcdf file.
>>>>
>>>> Not sure what I'm doing wrong here, printing the ncdf_out at the end
>>>> shows
>>>> only 1 item in the time dimension. Doesn't matter if I place the
>>>> following
>>>> code in the do loop or after it. I don't know how to associate the grid
>>>> with the time array and I couldn't understand how this should be done.
>>>> Also, the model I'm working with needs data to be dimensioned as:
>>>> y,x,time.
>>>>
>>>> Variable: ncdf_out
>>>> Type: file
>>>> filename: test
>>>> path: ./test.nc
>>>> file global attributes:
>>>> title : Temperature Data: degC
>>>> source : NLDAS Forcing2 - hourly
>>>> creation_date : Fri Mar 14 14:01:28 MDT 2014
>>>> dimensions:
>>>> lat = 19
>>>> lon = 17
>>>> time = 1
>>>> variables:
>>>> float lat ( lat )
>>>>
>>>> float lon ( lon )
>>>>
>>>> float time ( time )
>>>>
>>>> float TEMP ( lat, lon, time )
>>>> _FillValue : 9.96921e+36
>>>>
>>>>
>>>>
>>>> Many thanks in advance for your help,
>>>> john
>>>>
>>>>
>>>>
>>>>
>>>> Here's my code:
>>>> ;**************************************************
>>>> begin
>>>> ncfiles=systemfunc("ls *.nc")
>>>> if(dimsizes(ncfiles) .gt. 0) then
>>>> system("rm *.nc") ;remove any pre-existing file
>>>> end if
>>>>
>>>> ;specify extent
>>>> latN=42.852
>>>> latS=40.500
>>>> lonL=-112.620
>>>> lonR=-110.480
>>>>
>>>> cmd="ls *.grb"
>>>> grbFiles=systemfunc(cmd)
>>>> nfiles=dimsizes(grbFiles)
>>>>
>>>> ;create time var
>>>> ;should be 'hours since 2010-10-01T00.00
>>>> ;time should be provided as floating numbers
>>>> start=str_split(grbFiles(0),".")
>>>> year=str_get_cols(start(1), 1, 4)
>>>> month=str_get_cols(start(1), 5, 6)
>>>> day=str_get_cols(start(1), 7, 8)
>>>> hour=str_get_cols(start(2),0,1)
>>>>
>>>> time = fspan(0,nfiles,1)
>>>> time!0 = "time"
>>>> time@units = "hours since "+year+"-"+month+"-"+day+"T00.00"
>>>> time&time = time
>>>> timec=tocharacter(time)
>>>> ;print(time)
>>>>
>>>>
>>>> ;create and define the output file
>>>> ncdf_out = addfile("./test.nc" ,"c")
>>>> ncdf_out@title="Temperature Data: degC"
>>>> ncdf_out@source="NLDAS Forcing2 - hourly"
>>>> ncdf_out@creation_date = systemfunc ("date")
>>>> setfileoption(ncdf_out,"DefineMode",True)
>>>>
>>>>
>>>> ;get dimensions and setup file
>>>> g=addfile(grbFiles(0),"r")
>>>> l=g->lat_110
>>>> lat=l({latS:latN})
>>>> lat!0="lat"
>>>> lon=g->lon_110({lonL:lonR})
>>>> lon!0="lon"
>>>> nlat=dimsizes(lat)
>>>> nlon=dimsizes(lon)
>>>>
>>>> dimNames = (/"lat","lon","time"/)
>>>> dimSizes = (/nlat,nlon,dimsizes(time)/)
>>>> dimUnlim = (/False,False,False/)
>>>> filedimdef(ncdf_out,dimNames,dimSizes,dimUnlim)
>>>>
>>>>
>>>> ;setup the variable
>>>> t=g->TMP_110_HTGL
>>>> tc = new((/nlat,nlon,dimsizes(time)/), float)
>>>> tc!0="lat"
>>>> tc!1="lon"
>>>> tc!2="time"
>>>> tc&time=time
>>>> tc&lat=lat
>>>> tc&lon=lon
>>>> tc@sub_center=t@sub_center
>>>> tc@center=t@center
>>>> tc@long_name=t@long_name
>>>> tc@units="degC"
>>>> tc@level_indicator=t@level_indicator
>>>> tc@gds_grid_type=t@gds_grid_type
>>>> tc@parameter_table_version=t@parameter_table_version
>>>> tc@parameter_number=t@parameter_number
>>>> tc@model=t@model
>>>> tc@level=t@level
>>>>
>>>> ;setup file variables
>>>> filevardef(ncdf_out, "lat" ,typeof(lat),getvardims(lat))
>>>>
>>>> filevardef(ncdf_out, "lon" ,typeof(lon),getvardims(lon))
>>>> filevardef(ncdf_out, "time" ,typeof(time),getvardims(time))
>>>>
>>>> filevardef(ncdf_out, "TEMP" ,typeof(tc) ,getvardims(tc))
>>>>
>>>>
>>>> ; loop through the files and process
>>>> do n = 0, nfiles-1
>>>> g=addfile(grbFiles(n),"r")
>>>>
>>>> ;do the conversions and add to ncdf_out
>>>> ; TEMPERATURE
>>>> t=g->TMP_110_HTGL
>>>> tc=t({latS:latN},{lonL:lonR})-273.15
>>>> ;tc(,,n)=t({latS:latN},{lonL:lonR})-273.15
>>>> ;this spot reserved for regridding to higher resolution
>>>> ;ncdf_out->TEMP=(/tc/)
>>>>
>>>>
>>>> tc!0="lat"
>>>> tc!1="lon"
>>>> tc!2="time"
>>>> tc&lat=lat
>>>> tc&lon=lon
>>>> tc&time=time
>>>> end do
>>>>
>>>>
>>>> ncdf_out->TEMP=(/tc/)
>>>> print(ncdf_out)
>>>> end
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Mon Mar 17 07:58:04 2014

This archive was generated by hypermail 2.1.8 : Fri Mar 21 2014 - 15:49:21 MDT