[ncl-talk] Trouble converting dat to netcdf
Stanley Edwin
sedwin at alaska.edu
Sat Mar 19 15:43:01 MDT 2016
70200.dat.gz
<https://drive.google.com/a/alaska.edu/file/d/0B5e9G7If5LpKNnpnYXc0a3NfODg/view?usp=drive_web>
I could not get this dat file to convert properly to NetCDF, files from NWS
Thank You
Stanley G. Edwin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160319/62990b2f/attachment.html
-------------- next part --------------
;----------------------------------------------------------------------
; This script reads the 7 columns of data off a 70261.dat file and
; writes 6 columns of the data and the time stamps to a NetCDF file.
;----------------------------------------------------------------------
;---Read DAT file as strings so we can parse them.
file_prefix = "70200"
dat_filename = file_prefix + ".dat"
lines = asciiread(dat_filename,-1,"string")
;---Get index values of lines that contain a "#"
sects = str_match_ind(lines,"#")
nsect = dimsizes(sects)
;
; Loop through each subsection and store the # of rows for
; each section in an array "nrows".
;
; Also parse the timestamp and save the date info to
; individual arrays.
;
maxLevel = 0
nrows = new(nsect,integer)
yyyy = new(nsect,integer)
mm = new(nsect,integer)
dd = new(nsect,integer)
hh = new(nsect,integer)
do i=0,nsect-1
tmp = str_split(lines(sects(i))," ")
;print("tmp(0) = " + tmp(0) + ", tmp(1) = " + tmp(1))
yyyy(i) = toint(str_get_cols(tmp(0),6,9))
mm(i) = toint(str_get_cols(tmp(0),10,11))
dd(i) = toint(str_get_cols(tmp(0),12,13))
hh(i) = toint(str_get_cols(tmp(0),14,15))
nrows(i) = toint(tmp(1))
end do
npts = sum(nrows)
print("Total # npts = " + npts)
maxLevel = max(nrows)
print("MaxLevel = " + maxLevel)
;---Construct a time array
mn = new(nsect,integer)
sc = new(nsect,double)
mn = 0
sc = 0
units = "hours since " + min(yyyy) + "-01-01 00:00:00"
time = cd_inv_calendar(yyyy,mm,dd,hh,mn,sc,units, 0)
printMinMax(time,0)
;
; Loop through each section again, and this time
; read the data into the data_array variable.
;
; Each row of data is formatted into lined-up columns.
; For example:
;
;21 95700B 135 -290B-9999 0 0
;21100100B 135 -289B-9999 120 10
;
; Note that it is possible for columns to run into each other.
;
; Column #1: 00-01
; Column #2: 02-07 (column 8 is a letter)
; Column #3: 09-13 (column 14 is a letter)
; Column #4: 15-19 (column 20 is a letter)
; Column #5: 21-25
; Column #6: 26-30
; Column #7: 31-35
;
;---Create empty arrays to hold all columns of data.
p = new((/nsect, maxLevel/),float)
h = new((/nsect, maxLevel/),float)
t = new((/nsect, maxLevel/),float)
dp = new((/nsect, maxLevel/),float)
wdir = new((/nsect, maxLevel/),float)
wspd = new((/nsect, maxLevel/),float)
;
; Loop through each section of rows and parse the
; seven columns of data.
;
n = 0 ; counter
do i=0,nsect-1
ibeg = sects(i)+1 ; sect(i) is the "#xxxxx" line
iend = sects(i)+nrows(i)
p(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 2, 7))
h(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 9, 13))
t(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 15, 19))
dp(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 21, 25))
wdir(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 26, 30))
wspd(i, 0:nrows(i)-1) = tofloat(str_get_cols(lines(ibeg:iend), 31, 35))
end do
;---Make sure -9999 and -8888 values are flagged as missing.
p = where(p.eq.-9999.or.p.eq.-8888,p at _FillValue,p)
h = where(h.eq.-9999.or.h.eq.-8888,h at _FillValue,h)
t = where(t.eq.-9999.or.t.eq.-8888,t at _FillValue,t)
dp = where(dp.eq.-9999.or.dp.eq.-8888,dp at _FillValue,dp)
wdir = where(wdir.eq.-9999.or.wdir.eq.-8888,wdir at _FillValue,wdir)
wspd = where(wspd.eq.-9999.or.wspd.eq.-8888,wspd at _FillValue,wspd)
;print(nrows(0))
;print(p(0, :))
;print(h(0, :))
;print(t(0, :))
;print(dp(0, :))
;print(wdir(0, :))
;print(wspd(0, :))
;---Attach some metadata
time!0 = "time"
time&time = time ; This makes "time" a coordinate array
printVarSummary(time)
levs = ispan(0, maxLevel-1, 1)
levs!0 = "levs"
levs&levs = levs ; This makes "levs" a coordinate array
printVarSummary(levs)
;print(levs)
nrows!0 = "time"
nrows&time = time
nrows at long_name = "Available Levels"
nrows at units = "None"
p!0 = "time"
p!1 = "levs"
p at long_name = "pressure"
p at units = "Pa"
h!0 = "time"
h!1 = "levs"
h at long_name = "height"
h at units = "m"
t!0 = "time"
t!1 = "levs"
t at long_name = "temperature"
t at units = "degC"
dp!0 = "time"
dp!1 = "levs"
dp at long_name = "dew point"
dp at units = "degC"
wdir!0 = "time"
wdir!1 = "levs"
wdir at long_name = "wind direction"
wdir at units = "deg"
wspd!0 = "time"
wspd!1 = "levs"
wspd at long_name = "wind speed"
wspd at units = "m/s"
;---Print out min/max of each column
printMinMax(p,0)
printMinMax(h,0)
printMinMax(t,0)
printMinMax(dp,0)
printMinMax(wdir,0)
printMinMax(wspd,0)
;----------------------------------------------------------------------
; Section to write data to a NetCDF file
;----------------------------------------------------------------------
;---Make sure the file doesn't already exist.
cdf_filename = "NWS"+file_prefix + ".nc"
system("/bin/rm -f " + cdf_filename)
fout = addfile(cdf_filename, "c")
;---Write global attributes to the file (optional)
fout at title = "NCL Simple Approach to netCDF Creation"
fout at Conventions = "None"
fout at creation_date = systemfunc ("date")
filedimdef(fout, "time", -1, True)
;---Write data values (and metadata) to NetCDF file.
fout->p = p
fout->h = h
fout->t = t
fout->dp = dp
fout->wdir = wdir
fout->wspd = wspd
fout->time = time
fout->avlevs = nrows
More information about the ncl-talk
mailing list