Re: How to read gpcp daily data in NCL?

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri, 25 Sep 2009 11:57:45 -0600

I also have an NCL code that does this. It was written some time ago.
I made some minor changes. Mateus read the header and data records via
fortran.
I skipped the header but derived the year and month from the file name.

I ran the code to produce a file [603MB].

ftp ftp.cgd.ucar.edu
anonymous
email
cd /ftp/pub/shea/GPCP
binary
get GPCP.199610-200906.nc
quit
=======================================================

The data are in "short" .... use short2flt to unpack

netcdf GPCP.199610-200906 {
dimensions:
        time = UNLIMITED ; // (4656 currently)
        lon = 360 ;
        lat = 180 ;
variables:
        double time(time) ;
                time:long_name = "time" ;
                time:calendar = "standard" ;
                time:units = "days since 1990-01-01 00:00:00" ;
        float lat(lat) ;
                lat:long_name = "latitude" ;
                lat:units = "degrees_north" ;
        float lon(lon) ;
                lon:long_name = "longitude" ;
                lon:units = "degrees_east" ;
        int date(time) ;
                date:units = "yyyymmddhh" ;
                date:long_name = "gregorian date" ;
        short PREC(time, lat, lon) ;
                PREC:_FillValue = 32767s ;
                PREC:units = "mm/day" ;
                PREC:long_name = "GPCP: daily precipitation" ;
                PREC:missing_value = 32767s ;
                PREC:add_offset = 0.f ;
                PREC:scale_factor = 0.1f ;
                PREC:vMin_user_specified = 0.f ;
                PREC:vMax_user_specified = 2000.f ;
                PREC:vRange = 2000.f ;

// global attributes:
                :Title = "GPCP ONE-DEGREE DAILY PRECIPITATION DATA SET" ;
                :Comment = "netCDF version of original binary file(s)" ;
                :GSFC = "http://precip.gsfc.nasa.gov/" ;
                :Source = "ftp://rsd.gsfc.nasa.gov/pub/1dd-v1.1/" ;
                :Convention = "CF-1.0" ;
                :Conversion = "NCL: http://www.ncl.ucar.edu/" ;
                :Reference = "\n",
                        "Huffman, G.J., R.F. Adler, M.M. Morrissey, S.
Curtis \n",
                        "R. Joyce, B. McGavock, and J. Susskind,
2001: \n",
                        "Global precipitation at one-degree daily
resolution from multi-satellite observations\n",
                        "J. Hydrometeor., 2, 36-50\n",
                        "" ;
                :creation_date = "Fri Sep 25 11:29:33 MDT 2009" ;

Mateus Teixeira wrote:
> Hi Sabeer,
>
> I made this some time ago. Firstly, I create a Fortran code to take
> off the header of gpcp daily binary files. In such time I couldn't do
> this directly with NCL. So, I read them (the new files without header)
> within NCL and then create the netCDF file. I have attached the
> fortran program and the ncl script for you take a look.
>
> Unfortunately, the comments are in portuguese, but if you have any
> doubt you can just ask.
>
> Best regards,
>
> Mateus
>
>
> 2009/9/25 Sabeerali(sebi) <sabeerl_at_gmail.com>:
>
>> Helloo
>>
>> I am having GPCP daily data from 1996 to 2008....I want to read
>> those data and convert to netcdf...Anyone has done it earlier?
>>
>> Thanks in advance,
>> Sabeer
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>>
>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================

; -------------------------------------------------------
; GPCP ONE-DEGREE DAILY PRECIPITATION DATA SET
;
;http://precip.gsfc.nasa.gov/
;
; -------------------------------------------------------
; read the documentation for updates and details
; -------------------------------------------------------
; Convert binary to netCDF
; -------------------------------------------------------

  DIRI = "/project/cas/shea/GPCP/DailyBinary/" ; input directory with binary files
  DIRO = "/scratch/shea/" ; output directory with netCDF

  PACK = True
  if (PACK) then
      optPack = True
      optPack_at_min_value = 0.
      optPack_at_max_value = 2000. ; max value thru 2009 = 1975 mm/day :-)
      optPack_at_scale_factor = 0.1
      optPack_at_add_offset = 0.0
  else
      optPack = False
  end if

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

function parse_gpcp_filename( fName:string )
; parse file name: return yyyy, mm
; gpcp_1dd_v1.1_p1d.200906
; 012345678901234567890123
; 1 2
begin
  info = new( 2, "integer", "No_FillValue")
  c = stringtochar( fName )
  info(0) = stringtoint( (/c(18:21)/) ) ; yyyy
  info(1) = stringtoint( (/c(22:23)/) ) ; mm
  return( info )
end
  
begin

; -------------------------------------------------------
; define lat/lon grid: GPCP goes from N=>S
; -------------------------------------------------------
  nlat = 180
  mlon = 360
  lon = lonGlobeFo(mlon,"lon","longitude","degrees_east")
  lat = latGlobeFo(nlat,"lat","latitude","degrees_north")
  lat = lat(::-1) ; make N=>S
  printVarSummary(lat)

; -------------------------------------------------------
; Define generic file attributes
; -------------------------------------------------------

   nline = inttochar(10) ; new line character
   fAtt = 0 ; attributes for netCDF file
   fAtt_at_creation_date = systemfunc("date")

   fAtt_at_Reference = nline + \
       "Huffman, G.J., R.F. Adler, M.M. Morrissey, S. Curtis " + nline + \
       "R. Joyce, B. McGavock, and J. Susskind, 2001: " + nline + \
       "Global precipitation at one-degree daily resolution from multi-satellite observations" + nline +\
       "J. Hydrometeor., 2, 36-50" + nline

   fAtt@Conversion = "NCL: http://www.ncl.ucar.edu/"
   fAtt_at_Convention = "CF-1.0"
   fAtt@Source = "ftp://rsd.gsfc.nasa.gov/pub/1dd-v1.1/"
   fAtt@GSFC = "http://precip.gsfc.nasa.gov/"
   fAtt_at_Comment = "netCDF version of original binary file(s)"
   fAtt_at_Title = "GPCP ONE-DEGREE DAILY PRECIPITATION DATA SET"

; -------------------------------------------------------
; names of all 1dd files in the directory "diri"
; -------------------------------------------------------
   diri = DIRI
   fili = systemfunc("cd "+diri+" ; ls gpcp_1dd_v1.1_p1d*")
   nfil = dimsizes( fili )
   print(fili)

   diro = DIRO
; -------------------------------------------------------
; Loop over all files and create individual netCDF.
; Use "ncrcat" if one file is desired/
; -------------------------------------------------------
  setfileoption("bin","ReadByteOrder","BigEndian")

  do nf=0,nfil-1
     info = parse_gpcp_filename( fili(nf) ) ; yyyy, mm
     yyyy = info(0) ; clarity
     mm = info(1)
     ntim = days_in_month(yyyy, mm )

          ; read header + all days in the month
     dumy = fbindirread(diri+fili(nf),0,360+ntim*360*180,"float")
          ; create netCDF
     prec = onedtond( dumy(360:), (/ntim,nlat,mlon/) ) ; skip header
     delete (dumy) ; size may change for next month

; ---------------------------------------------------------
; Construct primary data object
; ---------------------------------------------------------
     prec_at_long_name = "GPCP: daily precipitation"
     prec_at_units = "mm/day"
     prec@_FillValue = -99999.
     prec_at_missing_value = prec@_FillValue

     prec!0 = "time"
     prec!1 = "lat"
     prec!2 = "lon"

     prec&lat = lat
     prec&lon = lon

; ---------------------------------------------------------
; create time variables
; ---------------------------------------------------------
     yyyymm = yyyy*100 + mm
     days = ispan(1,ntim,1)
     hh = 12 ; middle of day (arbitrary)

     date = yyyymm*10000 + days*100 + hh
     date!0 = "time"
     date_at_long_name = "gregorian date"
     date_at_units = "yyyymmddhh"
     
     YYYY = conform(date, yyyy, -1) ; make scalar a vector to match "days"
     MM = conform(date, mm, -1)
     HH = conform(date, hh, -1)
     ZERO = conform(date, 0, -1)

     tunits = "days since 1990-01-01 00:00:00" ; arbitrary
     time = ut_inv_calendar(YYYY,MM,days,HH, ZERO , ZERO ,tunits, 0)
     time!0 = "time"
     time_at_long_name = "time"
     time_at_units = tunits

; ---------------------------------------------------------
; create netCDF
; ---------------------------------------------------------
     yyyy = info(0)
     mm = info(1)
     yyyymm = yyyy*100 + mm
    ;ncfile = diro+fili(nf)+".nc"
     ncfile = diro+"GPCP.1DD."+yyyymm+".nc"
     print (ncfile)
     system ("/bin/rm -f " + ncfile) ; remove an pre-file
   
     ncdf = addfile(ncfile,"c") ; "c"reate the netCDF file
   
     fileattdef( ncdf, fAtt )
   
     dimNames = (/ "time", "lon", "lat" /)
     dimSizes = (/ -1 , mlon, nlat /)
     dimUnlim = (/ True , False, False /)
     filedimdef( ncdf, dimNames, dimSizes, dimUnlim )
                                            ; Define 1D variables.
     filevardef ( ncdf, "time", typeof(time), getvardims(time) )
     filevarattdef( ncdf, "time", time )
   
     filevardef ( ncdf, "lat", typeof(lat), getvardims(lat) )
     filevarattdef( ncdf, "lat", lat )
   
     filevardef ( ncdf, "lon", typeof(lon), getvardims(lon) )
     filevarattdef( ncdf, "lon", lon )
                                            ; Define 1D variables.
     filevardef ( ncdf, "date", typeof(date), getvardims(date) )
     filevarattdef( ncdf, "date", date )
   
     if (PACK) then
         filevardef ( ncdf, "PREC", "short", getvardims(prec) )
         pShort = pack_values(prec, "short", optPack)
         delete(pShort_at_vMin_original_data) ; extraneous
         delete(pShort_at_vMax_original_data)
         filevarattdef( ncdf, "PREC", pShort )
     else
         filevardef ( ncdf, "PREC", typeof(prec), getvardims(prec) )
         filevarattdef( ncdf, "PREC", prec )
     end if
                                               ; Write variables.
     ncdf->time = (/time/)
     ncdf->lat = (/lat /)
     ncdf->lon = (/lon /)
     ncdf->date = (/date/)
     if (PACK) then
         ncdf->PREC = (/pShort/)
         delete(pShort)
     else
         ncdf->PREC = (/prec/)
     end if

     delete (prec) ; size may change for next month
     delete (time)
     delete (date)
     delete (days)
     delete(YYYY )
     delete(MM )
     delete(HH )
     delete(ZERO )
  end do
  
end

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
   nt = 0
   pltType = "x11"
   pltName = "GPCP_DAILY"

   dira = "/scratch/shea/"
   fila = "GPCP.199610-200906.nc"

   f = addfile(dira+fila, "r")
   if (getfilevartypes(f, "PREC") .eq."short") then
       prec= short2flt( f->PREC )
   else
       prec= f->PREC
   end if

   printVarSummary(prec)
   printMinMax(prec,True)

   date = f->date

   res = True ; plot mods desired
   res_at_gsnMaximize = True
   
   yyyymm = date/100
   pltName = pltName+"."+yyyymm(nt) ; add date to name

   wks = gsn_open_wks(pltType, pltName)
   gsn_define_colormap(wks,"BlAqGrYeOrReVi200")

   res_at_gsnMaximize = True ; make ps/eps/pdf large
   res_at_cnFillOn = True
   res_at_cnFillMode = "RasterFill" ; Raster Mode
   res_at_cnLinesOn = False ; Turn off contour lines
   res_at_cnLineLabelsOn = False ; Turn off contour lines
   res_at_cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels
   if (prec_at_units.eq."mm/day") then
       res_at_cnLevels = (/0.5,1,2,3,4,5,6,7,8,9,10,20,30,50,100/) ; mm/day
   end if
   res_at_cnFillColors= (/ 0,38,47,55,70,86,94,102,116,127,142,163,178,186,192,201/)

   res_at_gsnCenterString = "GPCP: "+yyyymm(nt)

;;;prec = where(prec.eq.0 ,prec@_FillValue, prec)

   plot = gsn_csm_contour_map_ce(wks,prec(nt,:,:),res) ; create plot
   printMinMax(prec(nt,:,:),True)
end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Sep 25 2009 - 11:57:45 MDT

This archive was generated by hypermail 2.2.0 : Sun Sep 27 2009 - 11:32:26 MDT