Validating rain along the trajectory.

[midhun]

Hi Barbara,
Thanks!
Now it is clear

[yaqiang]

Please look at below sample script:

Code: Select all

#--------------------------------------------------------        
# Author: Yaqiang Wang                                           
# Date: 2016-5-6                                               
# Purpose: Convert ARL data to netCDF data  
# Note: Sample                                                   
#-----------------------------------------------------------

infn = 'D:/Temp/arl/gdas1.jul09.w5'    #Input ARL file
outfn = 'D:/Temp/arl/pr_gdas1.jul09.w5.nc'    #Output netCDF file

#Open ARL data file
f = addfile(infn)
rain = f['TPP6']
tnum = rain.dimlen(0)
ynum = rain.dimlen(2)
xnum = rain.dimlen(3)
lat = rain.dimvalue(2)
lon = rain.dimvalue(3)

#New netCDF file
ncfile = addfile(outfn, 'c')

#Add dimensions
xdim = ncfile.adddim('x', xnum)
ydim = ncfile.adddim('y', ynum)
tdim = ncfile.adddim('time', tnum)

#Add global attributes
ncfile.addgroupattr('Conventions', 'CF-1.6')
ncfile.addgroupattr('history', 'source: ftp://arlftp.arlhq.noaa.gov/pub/archives/gdas1/gdas1.jul09.w5')

#Add coordinate variables
yvar = ncfile.addvar('lat', 'float', [ydim])    #Latitude
yvar.addattr('long_name', 'Latitude')
yvar.addattr('units', 'degrees_north')
yvar.addattr('axis', 'Y')

xvar = ncfile.addvar('lon', 'float', [xdim])    #Longitude
xvar.addattr('long_name', 'Longitude')
xvar.addattr('units', 'degrees_east')
xvar.addattr('axis', 'X')

tvar = ncfile.addvar('time', 'int', [tdim])    #Time
tvar.addattr('long_name', 'Time')
tvar.addattr('units', 'hours since 1900-01-01 00:00:0.0')
tvar.addattr('axis', 'T')

#Add data variables
prvar = ncfile.addvar('PR', 'float', [tdim, ydim, xdim])

#Create netCDF file
ncfile.create()

#Write variable data
ncfile.write(xvar, lon)
ncfile.write(yvar, lat)
sst = datetime.datetime(1900,1,1)
for t in range(0, tnum):
    st = f.gettime(t)
    print st.strftime('%Y-%m-%d %H:00')
    hours = (st - sst).total_seconds() // 3600
    origin = [t]
    ncfile.write(tvar, array([hours]), origin=origin)
    print 'Variable: ' + prvar.name 
    data = f['TPP6'][t,:,:,:]              
    origin = [t, 0, 0]
    shape = [1, ynum, xnum]
    data = data.reshape(shape)
    ncfile.write(prvar, data, origin=origin)

#Close netCDF file
ncfile.close()

print 'Finished!'

[yaqiang]

I forgot to add attribute to precipitation variable in above script. See revised one:

Code: Select all

#--------------------------------------------------------        
# Author: Yaqiang Wang                                           
# Date: 2016-5-6                                               
# Purpose: Convert ARL data to netCDF data  
# Note: Sample                                                   
#-----------------------------------------------------------

infn = 'D:/Temp/arl/gdas1.jul09.w5'    #Input ARL file
outfn = 'D:/Temp/arl/pr_gdas1.jul09.w5.nc'    #Output netCDF file

#Open ARL data file
f = addfile(infn)
rain = f['TPP6']
tnum = rain.dimlen(0)
ynum = rain.dimlen(2)
xnum = rain.dimlen(3)
lat = rain.dimvalue(2)
lon = rain.dimvalue(3)

#New netCDF file
ncfile = addfile(outfn, 'c')

#Add dimensions
xdim = ncfile.adddim('x', xnum)
ydim = ncfile.adddim('y', ynum)
tdim = ncfile.adddim('time', tnum)

#Add global attributes
ncfile.addgroupattr('Conventions', 'CF-1.6')
ncfile.addgroupattr('history', 'source: ftp://arlftp.arlhq.noaa.gov/pub/archives/gdas1/gdas1.jul09.w5')

#Add coordinate variables
yvar = ncfile.addvar('lat', 'float', [ydim])    #Latitude
yvar.addattr('long_name', 'Latitude')
yvar.addattr('units', 'degrees_north')
yvar.addattr('axis', 'Y')

xvar = ncfile.addvar('lon', 'float', [xdim])    #Longitude
xvar.addattr('long_name', 'Longitude')
xvar.addattr('units', 'degrees_east')
xvar.addattr('axis', 'X')

tvar = ncfile.addvar('time', 'int', [tdim])    #Time
tvar.addattr('long_name', 'Time')
tvar.addattr('units', 'hours since 1900-01-01 00:00:0.0')
tvar.addattr('axis', 'T')

#Add data variables
prvar = ncfile.addvar('PR', 'float', [tdim, ydim, xdim])
prvar.addattr('long_name', 'Accumulated precipitation (6 h accumulation)')
prvar.addattr('units','m')

#Create netCDF file
ncfile.create()

#Write variable data
ncfile.write(xvar, lon)
ncfile.write(yvar, lat)
sst = datetime.datetime(1900,1,1)
for t in range(0, tnum):
    st = f.gettime(t)
    print st.strftime('%Y-%m-%d %H:00')
    hours = (st - sst).total_seconds() // 3600
    origin = [t]
    ncfile.write(tvar, array([hours]), origin=origin)
    print 'Variable: ' + prvar.name 
    data = f['TPP6'][t,:,:,:]              
    origin = [t, 0, 0]
    shape = [1, ynum, xnum]
    data = data.reshape(shape)
    ncfile.write(prvar, data, origin=origin)

#Close netCDF file
ncfile.close()

print 'Finished!'

[midhun]

Dear Yaqiang,
Thank you for the help.
It is a great tool!

[barbara.stunder]

The latitude and longitude of each gridpoint are not directly in the GDAS1 files. Instead the grid specification parameters are in the ASCII part of the index (INDX) record for each time period. The grid specs are in Table 1 at http://www.ready.noaa.gov/gdas1.php.
In linux 'head -c152 gda1.mar16.w1' gives the 1st 152 characters in the 1st record of the file.

The date/time are the first 8 characters of the record.

[danny67]

Dear Yaqiang,

Thank you!
Great! It worked!

Please see my script below. I have defined the lat/lon/time according to the description given at http://ready.arl.noaa.gov/gdas1.php. Is there any way to get the lat/lon/time values directly from the input file? That's why I wanted to see here for some details.

Code: Select all

infn =  '/home/midhun/Downloads/gdas1.jun13.w1'
outfn = '/home/midhun/Downloads/pr_gdas1.jun13.w1.nc'
inf= addfile(infn)
rain=inf['TPP6'][:,:,:,:]

#New netCDF file
ncfile = addfile(outfn, 'c')
#Add dimensions
x=ncfile.adddim('lon',360)
y=ncfile.adddim('lat',181)
t=ncfile.adddim('time',56)
ncfile.addgroupattr('history', 'source: ftp://arlftp.arlhq.noaa.gov/pub/archives/gdas1/gdas1.jun13.w1')
#create variables
lat = ncfile.addvar('lat', 'float', [y])    #Latitude
lat.addattr('long_name', 'Latitude')
lat.addattr('units', 'degrees_north')
lat.addattr('axis', 'Y')

lon = ncfile.addvar('lon', 'float', [x])    #Longitude
lon.addattr('long_name', 'Longitude')
lon.addattr('units', 'degrees_east')
lon.addattr('axis', 'X')

time = ncfile.addvar('time', 'int', [t])
time.addattr('long_name', 'time')
time.addattr('units', 'hours since 1900-01-01 00:00:0.0')
time.addattr('axis', 'T')

pr = ncfile.addvar('PR', 'float', [t,y,x])    #
pr.addattr('long_name', 'Accumulated precipitation (6 h accumulation)')
pr.addattr('units','m') 

xx=arange(0,360,1)
yy=arange(-90,91,1)
tt=arange(0,168,3)

#write variables
ncfile.create()
ncfile.write(lat,yy,origin=[0])
ncfile.write(lon,xx,origin=[0])
ncfile.write(time,tt,origin=[0])
ncfile.write(pr,rain,origin=[0,0,0])

Hello, I am also trying to get the lat/lon values direct, did you figure out if that's possible? how would one do it?