Fill the missing days and read data into one array

From: Wen.J.Qu <wen.j.qu_at_nyahnyahspammersnyahnyah>
Date: Fri Sep 07 2012 - 16:24:14 MDT

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

Received on Fri Sep 7 16:24:23 2012

This archive was generated by hypermail 2.1.8 : Tue Sep 11 2012 - 15:30:42 MDT