Hi, Mary,
I am trying to construct code based on your script, to read the data of all the stations for all years into one array. I am firstly reading the data of one station for one year into array "data", then trying to array_append_record to expand the array to the data of one station for all years into array "data_new", thirdly trying to array_append_record to expand the array to the data of all the stations for all years into array "data_all".
In doing this, I also considering that there are some entire months or entires years missing for some stations, so I change the code a little bit to deal with that.
But when I run the code, it gets an error that "array_append_record: non-record dimensions must be the same size". I can not find where the problem is. Would you please help me with that? Thanks a lot in advance.
The code is attached.
All my best!
Shawn
Wen.J.Qu
2012-09-07
########################################################################################################
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;----------------------------------------------------------------------
; Function to create new yyyymmdd array based on existing
; yyyymmdd array that is missing some days.
;
; If the existing array goes from 19720423 to 19750813, then the
; new array will cover every day from 19720401 to 19750831.
;----------------------------------------------------------------------
function create_new_yyyymmdd(year,month)
local ntotal, ny, nm, ey, em, ict, nvals
begin
;
; yearmoday has several days missing. We want to generate a new
; yearmoday array with no gaps.
;
;---Get largest size we need for new yearmoda array.
nvals = dimsizes(year)
ntotal = (year(nvals-1)-year(0)+1)*366
yearmoda_new = new(ntotal,string)
ny = year(0)
nm = month(0)
ey = year(nvals-1)
em = month(nvals-1)
ict = 0
do while((ny.lt.year(nvals-1)).or.(ny.eq.ey.and.nm.le.em))
dim = days_in_month(ny,nm)
yearmoda_new(ict:ict+dim-1) = tostring((ny*10000)+(nm*100) + \
ispan(1,dim,1))
ict = ict + dim ; Increment yearmoda_new counter
if(nm.lt.12) then ; Increment year and month as necessary
nm = nm + 1
else
ny = ny+1
nm = 1
end if
end do
return(yearmoda_new(0:ict-1))
end
;----------------------------------------------------------------------
; Create new (larger) array given the size, and a list of indexes
; to place the new values.
;----------------------------------------------------------------------
function create_new_array(x,nsize,ii)
begin
if(isatt(x,"_FillValue")) then
x_new = new(nsize,typeof(x),x@_FillValue)
else
x_new = new(nsize,typeof(x))
end if
x_new(ii) = x
return(x_new)
end
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
diri = "/work/jwang7/quwj/parallel/data/gsod/" ; data directory
station_list = diri+"/station_list.txt" ; station_list contains station numbers of the station of interest
STA = asciiread(station_list,-1,"integer")
;---Find the appropriate ASCII file to read in
do sta_number_loop = 450000,600000,1 ; the station numbers of the station of interest are within this range
if(any(STA.eq.sta_number_loop)) then
sta_number = sta_number_loop
;Debug print
printVarSummary(sta_number)
print(sta_number)
do year_loop = 1973,2012,1
Year = year_loop
fili = systemfunc ("cd "+diri+sta_number+"/" + " ; ls *")
sta = toint( str_get_cols(fili, 0, 5) )
YEAR = toint( str_get_cols(fili,13,16) )
;Debug print
printVarSummary(Year)
print(Year)
printVarSummary(fili)
print(fili)
printVarSummary(sta)
print(sta)
printVarSummary(YEAR)
print(YEAR)
Index = ind(sta.eq.sta_number .and. YEAR.eq.Year)
if(ismissing(Index)) then
delete(fili)
delete(sta)
delete(YEAR)
delete(Index)
;Creat a new array for the station that year with missing values
nDayYear = 365
if (isleapyear(Year)) then
nDayYear = 366
end if
Yearmoda = yyyymmdd_time(Year,Year,"string")
delete(Year)
;Break yearmoda into individual year, month and day
year = toint(str_get_cols(Yearmoda,0,3))
month = toint(str_get_cols(Yearmoda,4,5))
day = toint(str_get_cols(Yearmoda,6,7))
year@calendar = "standard"
yearmoda_new = Yearmoda
stn_new = new(nDayYear,"integer",999999)
stn_new = sta_number
wban_new = new(nDayYear,"integer",99999)
temp_new = new(nDayYear,"float",9999.9)
dewp_new = new(nDayYear,"float",9999.9)
rh_new = new(nDayYear,"float",99.99)
slp_new = new(nDayYear,"float",9999.9)
stp_new = new(nDayYear,"float",9999.9)
visib_new = new(nDayYear,"float",999.9)
wdsp_new = new(nDayYear,"float",999.9)
maxspd_new = new(nDayYear,"float",999.9)
gust_new = new(nDayYear,"float",999.9)
maxtemp_new = new(nDayYear,"float",9999.9)
mintemp_new = new(nDayYear,"float",9999.9)
prcp_new = new(nDayYear,"float",99.99)
sndp_new = new(nDayYear,"float",999.9)
frshtt_new = new(nDayYear,"string")
fog_new = toint(str_get_cols(frshtt_new,0,0))
rain_new = toint(str_get_cols(frshtt_new,1,1))
snow_new = toint(str_get_cols(frshtt_new,2,2))
hail_new = toint(str_get_cols(frshtt_new,3,3))
thunder_new = toint(str_get_cols(frshtt_new,4,4))
tornado_new = toint(str_get_cols(frshtt_new,5,5))
yearmoda_new_int = toint(yearmoda_new)
frshtt_new_int = toint(frshtt_new)
;Construct a new array include all the varibles
data = (/stn_new,wban_new,yearmoda_new_int,temp_new,\
dewp_new,rh_new,slp_new,stp_new,visib_new,\
wdsp_new,maxspd_new,gust_new,maxtemp_new,\
mintemp_new,prcp_new,sndp_new,frshtt_new_int,\
fog_new,rain_new,snow_new,hail_new,thunder_new,\
tornado_new/)
;---Release varible
delete(year)
delete(month)
delete(day)
delete(Yearmoda)
delete(yearmoda_new)
delete(yearmoda_new_int)
delete(stn_new)
delete(wban_new)
delete(temp_new)
delete(dewp_new)
delete(rh_new)
delete(slp_new)
delete(stp_new)
delete(visib_new)
delete(wdsp_new)
delete(maxspd_new)
delete(gust_new)
delete(maxtemp_new)
delete(mintemp_new)
delete(prcp_new)
delete(sndp_new)
delete(frshtt_new)
delete(frshtt_new_int)
delete(fog_new)
delete(rain_new)
delete(snow_new)
delete(hail_new)
delete(thunder_new)
delete(tornado_new)
else
;Debug print
print (Index)
;Read the file's data
;Debug print
printVarSummary(fili(Index))
print(fili(Index))
filename = diri+sta_number+"/"+fili(Index)
;Release varibles
delete(fili)
delete(sta)
delete(YEAR)
delete(Index)
;Creat a new wholes array for the station that year
nDayYear = 365
if (isleapyear(Year)) then
nDayYear = 366
end if
Yearmoda = yyyymmdd_time(Year,Year,"string")
delete(Year)
;Break yearmoda into individual year, month and day
year = toint(str_get_cols(Yearmoda,0,3))
month = toint(str_get_cols(Yearmoda,4,5))
day = toint(str_get_cols(Yearmoda,6,7))
year@calendar = "standard"
;---Read data file as array of strings
lines = asciiread(filename,-1,"string")
nlines = dimsizes(lines)
;
; Parse strings to get values. Don't use header line [line(0)]
; Note that fields 5, 7, 9, 11, 13, and 15 are skipped.
; They don't appear to have headers associated with them.
;
stn = toint(str_get_field(lines(1:),1," "))
wban = toint(str_get_field(lines(1:),2," "))
yearmoda = str_get_field(lines(1:),3," ")
temp = tofloat(str_get_field(lines(1:),4," "))
dewp = tofloat(str_get_field(lines(1:),6," "))
slp = tofloat(str_get_field(lines(1:),8," "))
stp = tofloat(str_get_field(lines(1:),10," "))
visib = tofloat(str_get_field(lines(1:),12," "))
wdsp = tofloat(str_get_field(lines(1:),14," "))
maxspd = tofloat(str_get_field(lines(1:),16," "))
gust = tofloat(str_get_field(lines(1:),17," "))
maxtemp = tofloat(str_get_field(lines(1:),18," "))
mintemp = tofloat(str_get_field(lines(1:),19," "))
prcp = tofloat(str_get_field(lines(1:),20," "))
sndp = tofloat(str_get_field(lines(1:),21," "))
frshtt = str_get_field(lines(1:),22," ")
; Convert temperature from Fahrenheit to Celsius.
temp = (temp-32)*5/9
; Convert dew point temperature from Fahrenheit to Celsius.
dewp = (dewp-32)*5/9
; Calculate relative humidity (%).
rh = 100*(((112-0.1*temp+dewp)/(112+0.9*temp))^8)
; convert visib from mile to kilometer.
visib = visib*1.60934398
; convert wdsp from knot to m/s.
wdsp = wdsp*0.51444
; convert maxspd from knot to m/s.
maxspd = maxspd*0.51444
; convert gust from knot to m/s.
gust = gust*0.51444
; Convert maximum temperature from Fahrenheit to Celsius.
maxtemp = (maxtemp-32)*5/9
; Convert minimum temperature from Fahrenheit to Celsius.
mintemp = (mintemp-32)*5/9
; convert prcp from inch to millmeter.
prcp = prcp*25.4
; convert sndp from inch to millmeter.
sndp = sndp*25.4
nvals = dimsizes(stn)
;---Set some missing values where appropriate
stn@_FillValue = 999999
wban@_FillValue = 99999
temp@_FillValue = 9999.9
dewp@_FillValue = 9999.9
rh@_FillValue = 99.99
slp@_FillValue = 9999.9
stp@_FillValue = 9999.9
visib@_FillValue = 999.9
wdsp@_FillValue = 999.9
maxspd@_FillValue = 999.9
gust@_FillValue = 999.9
maxtemp@_FillValue = 9999.9
mintemp@_FillValue = 9999.9
prcp@_FillValue = 99.99
sndp@_FillValue = 999.9
;---Debug prints
print("stn min/max = " + min(stn) + "/" + max(stn))
print("wban min/max = " + min(wban) + "/" + max(wban))
print("temp min/max = " + min(temp) + "/" + max(temp))
print("dewp min/max = " + min(dewp) + "/" + max(dewp))
print("rh min/max = " + min(rh) + "/" + max(rh))
print("slp min/max = " + min(slp) + "/" + max(slp))
print("stp min/max = " + min(stp) + "/" + max(stp))
print("visib min/max = " + min(visib) + "/" + max(visib))
print("wdsp min/max = " + min(wdsp) + "/" + max(wdsp))
print("maxspd min/max = " + min(maxspd) + "/" + max(maxspd))
print("gust min/max = " + min(gust) + "/" + max(gust))
print("maxtemp min/max = " + min(maxtemp) + "/" + max(maxtemp))
print("mintemp min/max = " + min(mintemp) + "/" + max(mintemp))
print("prcp min/max = " + min(prcp) + "/" + max(prcp))
print("sndp min/max = " + min(sndp) + "/" + max(sndp))
printVarSummary(frshtt)
;---Create new array based on year and month arrays
yearmoda_new = create_new_yyyymmdd(year,month)
ntime = dimsizes(yearmoda_new)
print("--------------------------------------------------")
print("# of days in new arrays will be = " + ntime)
print("--------------------------------------------------")
;--Get indices where date value is not missing.
ii = new(nvals,integer)
do i=0,nvals-1
ii(i) = ind(yearmoda(i).eq.yearmoda_new)
end do
;---Use indexes to place real values in location of new arrays
stn_new = create_new_array(stn,ntime,ii)
wban_new = create_new_array(wban,ntime,ii)
temp_new = create_new_array(temp,ntime,ii)
dewp_new = create_new_array(dewp,ntime,ii)
rh_new = create_new_array(rh,ntime,ii)
slp_new = create_new_array(slp,ntime,ii)
stp_new = create_new_array(stp,ntime,ii)
visib_new = create_new_array(visib,ntime,ii)
wdsp_new = create_new_array(wdsp,ntime,ii)
maxspd_new = create_new_array(maxspd,ntime,ii)
gust_new = create_new_array(gust,ntime,ii)
maxtemp_new = create_new_array(maxtemp,ntime,ii)
mintemp_new = create_new_array(mintemp,ntime,ii)
prcp_new = create_new_array(prcp,ntime,ii)
sndp_new = create_new_array(sndp,ntime,ii)
frshtt_new = create_new_array(frshtt,ntime,ii)
fog_new = toint(str_get_cols(frshtt_new,0,0))
rain_new = toint(str_get_cols(frshtt_new,1,1))
snow_new = toint(str_get_cols(frshtt_new,2,2))
hail_new = toint(str_get_cols(frshtt_new,3,3))
thunder_new = toint(str_get_cols(frshtt_new,4,4))
tornado_new = toint(str_get_cols(frshtt_new,5,5))
printVarSummary(yearmoda_new)
yearmoda_new_int = toint(yearmoda_new)
printVarSummary(yearmoda_new_int)
frshtt_new_int = toint(frshtt_new)
;Construct a new array include all the varibles
;Debug print
printVarSummary(stn_new)
printVarSummary(wban_new)
printVarSummary(yearmoda_new_int)
printVarSummary(temp_new)
printVarSummary(dewp_new)
printVarSummary(rh_new)
printVarSummary(slp_new)
printVarSummary(stp_new)
printVarSummary(visib_new)
printVarSummary(wdsp_new)
printVarSummary(maxspd_new)
printVarSummary(gust_new)
printVarSummary(maxtemp_new)
printVarSummary(mintemp_new)
printVarSummary(prcp_new)
printVarSummary(sndp_new)
printVarSummary(frshtt_new_int)
printVarSummary(fog_new)
printVarSummary(rain_new)
printVarSummary(snow_new)
printVarSummary(hail_new)
printVarSummary(thunder_new)
printVarSummary(tornado_new)
data = (/stn_new,wban_new,yearmoda_new_int,temp_new,\
dewp_new,rh_new,slp_new,stp_new,visib_new,\
wdsp_new,maxspd_new,gust_new,maxtemp_new,\
mintemp_new,prcp_new,sndp_new,frshtt_new_int,\
fog_new,rain_new,snow_new,hail_new,thunder_new,\
tornado_new/)
;---Debug prints
print("stn_new min/max = " + min(stn_new) + "/" + max(stn_new))
print("wban_new min/max = " + min(wban_new) + "/" + max(wban_new))
print("temp_new min/max = " + min(temp_new) + "/" + max(temp_new))
print("dewp_new min/max = " + min(dewp_new) + "/" + max(dewp_new))
print("rh_new min/max = " + min(rh_new) + "/" + max(rh_new))
print("slp_new min/max = " + min(slp_new) + "/" + max(slp_new))
print("stp_new min/max = " + min(stp_new) + "/" + max(stp_new))
print("visib_new min/max = " + min(visib_new) + "/" + max(visib_new))
print("wdsp_new min/max = " + min(wdsp_new) + "/" + max(wdsp_new))
print("maxspd_new min/max = " + min(maxspd_new) + "/" + max(maxspd_new))
print("gust_new min/max = " + min(gust_new) + "/" + max(gust_new))
print("maxtemp_new min/max = " + min(maxtemp_new) + "/" + max(maxtemp_new))
print("mintemp_new min/max = " + min(mintemp_new) + "/" + max(mintemp_new))
print("prcp_new min/max = " + min(prcp_new) + "/" + max(prcp_new))
print("sndp_new min/max = " + min(sndp_new) + "/" + max(sndp_new))
printVarSummary(frshtt_new)
printVarSummary(data)
;---Release varible
delete(i)
delete(ii)
delete(year)
delete(month)
delete(day)
delete(lines)
delete(nlines)
delete(stn)
delete(wban)
delete(yearmoda)
delete(Yearmoda)
delete(temp)
delete(dewp)
delete(rh)
delete(slp)
delete(stp)
delete(visib)
delete(wdsp)
delete(maxspd)
delete(gust)
delete(maxtemp)
delete(mintemp)
delete(prcp)
delete(sndp)
delete(frshtt)
delete(yearmoda_new)
delete(yearmoda_new_int)
delete(stn_new)
delete(wban_new)
delete(temp_new)
delete(dewp_new)
delete(rh_new)
delete(slp_new)
delete(stp_new)
delete(visib_new)
delete(wdsp_new)
delete(maxspd_new)
delete(gust_new)
delete(maxtemp_new)
delete(mintemp_new)
delete(prcp_new)
delete(sndp_new)
delete(frshtt_new)
delete(frshtt_new_int)
delete(fog_new)
delete(rain_new)
delete(snow_new)
delete(hail_new)
delete(thunder_new)
delete(tornado_new)
end if
;Construct a new array for one station all the years
data_1 = data
;Release varible
delete(data)
if(year_loop.le.1973) then
data_new = data_1
;Release varible
delete(data_1)
else
data_2 = data_new
;Release varible "data_new"
delete(data_new)
data_new=array_append_record(data_1,data_2,0)
;Release varibles
delete(data_1)
delete(data_2)
end if
;Debug print
printVarSummary(data_new)
end do
;Construct an entire array for all the stations all the years
data_all_1 = data_new
;Release varible
delete(data_new)
if(sta_number_loop.le.450070) then
data_all = data_all_1
;Release varible
delete(data_all_1)
else
data_all_2 = data_all
;Release varible "data_all"
delete(data_all)
data_all=array_append_record(data_all_1,data_all_2,0)
;Release varibles
delete(data_all_1)
delete(data_all_2)
end if
;Debug print
printVarSummary(data_all)
end if
end do
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Tue Sep 11 2012 - 15:30:42 MDT