Here is some pretty efficient code I recently used to convert ECMWF
ERA40 data from grib to netcdf.
It is a modified version of some code originally sent to us by Dennis
Shea.
Hope this helps.
Ryan Sriver
;***************************************************
; era40 full variable conversion to netcdf
;***************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
function rdVar3D (f:file, vName:string, time:integer)
begin
var = f->$vName$ ; read variable
printVarSummary(var)
var!0 = "time" ; assign 'conventional' dimension names
var!1 = "lat"
var!2 = "lon"
; following not needed for current application
;var&time = time ; assign integer time coord variable
return(var)
end
begin
;***************************************************
; get all file names; then loop over each file
;***************************************************
diri = "./surface_state_4xdaily/" ; input directory
fils = systemfunc("cd "+diri+"; ls e4oper.an.sfc.200*")
nfils = dimsizes(fils)
print (fils)
OS = systemfunc ("uname")
diro = "./"
;***************************************************
; specify desired GRIB variables
;***************************************************
;system("banner HI")
do nf=0,nfils-1
print ("===========> nf="+nf+" <========================")
wcStrt = systemfunc("date") ; wall clock start
;***************************************************
; open grib file and get variable names
;***************************************************
grib_in = addfile(diri+fils(nf)+".grb","r")
names = getfilevarnames(grib_in) ; get ALL variable names
wallClockElapseTime(wcStrt, "Read/open file file: "+fils(nf), 0)
ncl_names = names
lat = grib_in->g4_lat_1 ; latitude
lat!0 = "lat" ; rename NCL's GRIB dimension
names
lon = grib_in->g4_lon_2 ; longitude
lon!0 = "lon" ; to 'conventional' dimension
stime = grib_in->initial_time0
time = grib_stime2itime (stime) ; string to integer
(contributed.ncl)
if (nf.eq.0) then
print (names)
print (stime+" "+time)
end if
ntim = dimsizes(time)
nlat = dimsizes(lat)
mlon = dimsizes(lon)
;***************************************************
; open output netcdf file
; create parameters and output file
;***********************************************
system ("rm -f "+diro+fils(nf)+".nc") ; remove pre-exit file
(if any)
netcdf_out = addfile(diro+fils(nf)+".nc","c") ; "c"reate the
netCDF file
;***********************************************
; assign file attributes
;***********************************************
fAtt = True
fAtt@title = "GRIB-to-netCDF: processed by
huberm@purdue.edu"
fAtt@source_file = fils(nf)
fAtt@Conventions = "None"
fAtt@System = OS
fAtt@creation_date = systemfunc ("date")
fileattdef( netcdf_out, fAtt )
;***********************************************
; predefine coordinate information
;***********************************************
dimNames = (/"time", "lat", "lon" /)
dimSizes = (/ ntim , nlat, mlon/)
dimUnlim = (/ True , False, False/)
filedimdef(netcdf_out, dimNames , dimSizes, dimUnlim )
filevardef (netcdf_out, "time" , typeof(time), "time" )
filevarattdef(netcdf_out, "time", time)
filevardef (netcdf_out, "lat", typeof(lat), "lat")
filevarattdef(netcdf_out, "lat", lat)
filevardef (netcdf_out, "lon", typeof(lon), "lon")
filevarattdef(netcdf_out, "lon", lon)
;***********************************************
; predefine variable sizes
; (unfortunately) must read variables from GRIB file
; needed for filevarattdef
;***********************************************
do i=3, dimsizes(ncl_names)-1
TMP = rdVar3D (grib_in, names(i) ,time)
filevardef(netcdf_out, names(i) \
, typeof(TMP), (/"time","lat","lon"/) )
filevarattdef(netcdf_out, names(i) , TMP)
end do
wallClockElapseTime(wcStrt, "Defintion phase took ", 0)
;***************************************************
; write data values to predefined locations
;***************************************************
wcStrt = systemfunc("date") ; wall clock start
netcdf_out->time = (/ time /)
netcdf_out->lat = (/ lat /)
netcdf_out->lon = (/ lon /)
do i=3, dimsizes(names)-1
TMP = rdVar3D (grib_in, names(i) ,time)
netcdf_out->$names(i)$ = (/ TMP /)
end do
wallClockElapseTime(wcStrt, "Write phase took ", 0)
delete (stime) ; may change next "nf" iteration
delete (time)
delete (names)
delete (lat)
delete (lon)
delete (TMP)
end do
end
On Feb 16, 2005, at 1:25 PM, David Brandon wrote:
> Does anyone have a good, GRIB to netcdf decoder. We now are getting
>
> Reanalyzed climate data in GRIB but want it in another format.
>
>
>
> Dave Brandon
>
> Hydrologist in Charge
>
> Colorado Basin River Forecast Center
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk@ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
Ryan L. Sriver
Department of Earth and Atmospheric Sciences
Purdue University
550 Stadium Mall Drive
West Lafayette, IN 47907-2051
Office: Civil Engineering Building Rm 4162
Phone: (765) 494-0652
e-mail: rsriver@purdue.edu
homepage: http://www.physics.purdue.edu/~rsriver
_______________________________________________
ncl-talk mailing list
ncl-talk@ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2b29 : Wed Feb 16 2005 - 14:20:29 MST