Re: GPCP 1dd binary read

Bonjour Prince,

I really wish the GPCP project would create netCDF. I don't like
the combined ascii header and binary data. Also, users
would not have to read external documentation. :-(

[1] "fbindirread" acts on files where *all* records have the *same* length.
      The GPCP daily grids are preceded by a 1440 character header.
       The daily grids follow. Hence, specifying the record 1 as you did
       will not work. The trick is the read the whole file and
       ignore the first 360 elements which equates to 1440 characters.
       [Note: 1440=4*360].

       The documentation and sample fortran codes mention this.

[2] I did have a script that takes the daily files and makes netCDF.
      Note: the following in the binary read
             360 + ntim*nlat*mlon
      The 360 preceding the plus sign is the number of "float" values
      that corresponds to 1440 characters. The GPCP
      creators designed the header of the is size specifically for this
purpose.

      I will send the full script to you offline. However, for ncl-talk
      the relevent statements are:

; -------------------------------------------------------
; read the documentation for details
; http://www1.ncdc.noaa.gov/pub/data/gpcp/1dd/doc/1DD_doc
; -------------------------------------------------------

function parse_gpcp_1dd_yyyymm( fName:string )
; parse file name: return yyyy, mm and number of days in month
begin
  date = new( 3, "integer", "No_FillValue")
  c = stringtochar( fName )
  date(0) = stringtoint( (/c(13:16)/) ) ; yyyy
  date(1) = stringtoint( (/c(17:18)/) ) ; mm
  date(2) = days_in_month ( date(0), date(1) ) ; number days in mm
  return( date )
end

nlat = 180
  mlon = 360
  diri = "/project/cas/shea/GPCP/"
  fili = systemfunc("cd "+diri+" ; ls gpcp_1dd_p1d.*")
  nfil = dimsizes( fili )

  setfileoption("bin","ReadByteOrder","BigEndian")

  do nf=0,nfils-1
     date = parse_gpcp_1dd_yyyymm( fili(nf) )
     ntim = date(2) ;
number of days in month

          ; read header + all days in the month
     dumy = fbindirread(firi_fili(nf),0,360+ntim*nlat*mlon,"float")
     prec = onedtond( dumy(360:), (/ntim,nlat,mlon/) )
     delete (dumy) ; size may change for next month
     :
     :
     :
     delete(prec) ; size may change for next month
  end do

Prince K. XAVIER wrote:
> Hi all,
>
> I am trying to read the GPCP 1DD data at
> http://www1.ncdc.noaa.gov/pub/data/gpcp/1dd/data/
> using the options:
>
> setfileoption("bin","ReadByteOrder","BigEndian")
> prec = fbindirread("gpcp_1dd_p1d.199707",1,(/180,360/),"float")
>
> But I find get the first record empty. I think I am doing something
> wrong with the record length. If someone has read this data using NCL,
> could you help to set the correct file options?
>
> thanks in advance
> Prince
>

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 26 2007 - 12:10:15 MST