Re: how to convert TRMM hdf data to nc files

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri, 12 Jan 2007 15:12:03 -0700 (MST)

Hello,

I was just about to delete the original email.

Then, I thought I'd take another look at this.
I have worked with TRMM 3B43V6 data.

I took a quick look at the 3B42 URL and downloaded
the following:

3B42.060101.0.6.HDF
3B42.060101.12.6.HDF
3B42.060101.15.6.HDF
3B42.060101.18.6.HDF
3B42.060101.21.6.HDF
3B42.060101.3.6.HDF
3B42.060101.6.6.HDF
3B42.060101.9.6.HDF

Note: "ls 3B42*HDF" will yield this order.
      This is *NOT* chronological. They *should* have created.

      3B42.060101.00.6.HDF
      3B42.060101.03.6.HDF
      3B42.060101.06.6.HDF
      3B42.060101.09.6.HDF
      3B42.060101.12.6.HDF
      3B42.060101.15.6.HDF
      3B42.060101.18.6.HDF
      3B42.060101.21.6.HDF

      Oh well!
      ==================================
I took an existing script and mangled it a bit
[TRMM.3B42.hdf2nc.ncl]

Thic create one netCDF for each input HDF.
I created an unlimited time dimension. The result is

-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010100.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010103.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010106.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010109.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010112.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010115.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010118.nc
-rw-r--r-- 1 shea cgdcas 2312532 Jan 12 14:58 3B42.2006010121.nc

these files can be read into ncl and then use addfile_getVar

or

%> ncrcat 3B42.2006*nc 3B42.20060101.nc
%> ncdump -v time,date,datesec 3B42.20060101.nc

This yielded

netcdf 3B42.20060101 {
dimensions:
        time = UNLIMITED ; // (8 currently)
        lat = 400 ;
        lon = 1440 ;
variables:
        double time(time) ;
                time:calendar = "gregorian" ;
                time:units = "hours since 1998-01-01 00:00:0.0" ;
        float lat(lat) ;
                lat:units = "degrees_north" ;
                lat:long_name = "latitude" ;
        float lon(lon) ;
                lon:units = "degrees_east" ;
                lon:long_name = "longitude" ;
        int date(time) ;
                date:units = "yyyymmdd" ;
                date:long_name = "current date" ;
        int datesec(time) ;
                datesec:units = "seconds" ;
                datesec:long_name = "current seconds of current date" ;
        float precip(time, lat, lon) ;
                precip:missing_value = -9999.f ;
                precip:_FillValue = -9999.f ;
                precip:avg_period = "0000-00-00 01:00:00" ;
                precip:delta_t = "0000-00-00 03:00:00" ;
                precip:units = "mm" ;
                precip:hdf_name = "precipitation" ;

// global attributes:
                :creation_date = "Fri Jan 12 14:58:25 MST 2007" ;
                :Conventions = "None" ;
                :source_html = "\n",
                        
"http://disc.sci.gsfc.nasa.gov/data/datapool/TRMM/01_Data_Products/02_Gridde
d/06_3-hour_Gpi_Cal_3B_42/index.html\n",
                        "" ;
                :title = "TRMM: 3B42" ;
                :history = "Fri Jan 12 14:59:09 2007: ncrcat 3B42.2006010100.nc
3B42.2006010103.nc 3
B42.2006010106.nc 3B42.2006010109.nc 3B42.2006010112.nc 3B42.2006010115.nc
3B42.2006010118.nc 3B42.2
006010121.nc 3B42.20060101.nc" ;
data:

 time = 70128, 70131, 70134, 70137, 70140, 70143, 70146, 70149 ;

 date = 20060101, 20060101, 20060101, 20060101, 20060101, 20060101, 20060101,
    20060101 ;

 datesec = 0, 10800, 21600, 32400, 43200, 54000, 64800, 75600 ;
}

good luck

function parse_3B42_date (fHDF[1]:string)
; parse a string with the following structure: 3B42.060101.0.6.HDF
; 3B42.060101.0.6.HDF ==> 3B42.yymmdd.h.6.HDF
; 3B42.060101.12.6.HDF ==> 3B42.yymmdd.hh.6.HDF
; return with a better file name==> 3B42.yyyymmddhh.6.nc
local perc, per, fhdfc
begin
  perc = stringtochar(".") ; dimsizes(chr)=2
  per = perc(0) ; period

  fhdfc = stringtochar( fHDF )
  iper = ind(fhdfc.eq.per)
                       
  yymmdd = stringtointeger( chartostring(fhdfc(iper(0)+1:iper(1)-1)) )
  hh = stringtointeger( chartostring(fhdfc(iper(1)+1:iper(2)-1)) )
  yy = yymmdd/10000
  if (yy.ge.90) then
      yyyy = 1900 + yy
  else
      yyyy = 2000 + yy
  end if
  mmdd = yymmdd - yy*10000
  mm = mmdd/100
  dd = mmdd - mm*100

  date = new ( 4, "integer", "No_FillValue")
  date(0) = yyyy
  date(1) = mm
  date(2) = dd
  date(3) = hh

  return(date)
end

begin ; begin not required

;*****************************************************
; Create TRMM coordinate variables
;*****************************************************
   nlat = 400
   mlon = 1440
   lat = ispan(0,nlat-1,1)*0.25 - 49.875
   lon = ispan(0,mlon-1,1)*0.25 - 179.875
  
   lat!0 = "lat"
   lat&lat = lat
   lat_at_long_name = "latitude"
   lat_at_units = "degrees_north"
 
   lon!0 = "lon"
   lon&lon = lon
   lon_at_long_name = "longitude"
   lon_at_units = "degrees_east"
                                      ; miscellaneous
   tunits = "hours since 1998-01-01 00:00:0.0"
   minute = 0
   second = 0.d0
   filAtt = True
   nline = inttochar(10)

;*****************************************************
; get all 3B42 HDF file names
; As returned by "ls", these will not be in chronological order
; files were poorly named
;*****************************************************
   diri = "./"
   fili = systemfunc("cd "+diri+" ; ls 3B42.*HDF")
   print(fili) ; these may not be in chronological oerder
   nfil = dimsizes( fili )
;*****************************************************
; loop over each file and create netCDF
;*****************************************************
   do n=0,nfil-1

      f = addfile ( fili(n), "r")

      P = f->precipitation
      P!0 = "time" ; rename to
      P!1 = "lon" ; match model names
      P!2 = "lat" ; convenience only

      p = P(time|:,lat|:,lon|:) ; reorder to match model

                                     ; ADD META DATA
      p_at_units = "mm" ; ?
      p_at_delta_t = "0000-00-00 03:00:00"
      p_at_avg_period = "0000-00-00 01:00:00"
      p@_FillValue = -9999.
      p_at_missing_value = p@_FillValue

                                     ; ADD COORDINATE VARIABLES
      p&lat = lat
      p&lon = lon

      ymdh = parse_3B42_date ( fili(n) )
      time = ut_inv_calendar(ymdh(0),ymdh(1),ymdh(2),ymdh(3) \
                              ,minute,second,tunits, 0)
      time_at_calendar = "gregorian"
      time!0 = "time"
      p&time = time
     ;printVarSummary(p)

                                     ; ADD AUXILARY TIME VARIABLES
      date = ymdh(0)*10000+ ymdh(1)*100 + ymdh(2)
      date!0 = "time"
      date_at_long_name = "current date"
      date_at_units = "yyyymmdd"

      datesec = ymdh(3)*3600
      datesec!0 = "time"
      datesec_at_long_name = "current seconds of current date"
      datesec_at_units = "seconds"

      yyyymmddhh = date*100 + ymdh(3)

    ;===================================================================
    ; WRITE OUTPUT netCDF: Use efficient method [ie: classic method]
    ;===================================================================
      ncDir = "./"
      ncName = "3B42."+yyyymmddhh+".nc"
      system ("'rm' -f "+ncDir+ncName)
      fnc = addfile (ncDir+ncName, "c")

    ;===================================================================
    ; explicitly declare file definition mode. Improve efficiency.
    ;===================================================================
        setfileoption(fnc,"DefineMode",True)

    ;===================================================================
    ; create global attributes of the file
    ;===================================================================
      filAtt_at_title = "TRMM: 3B42"
      filAtt_at_source_html = nline + \
"http://disc.sci.gsfc.nasa.gov/data/datapool/TRMM/01_Data_Products/02_Gridded/06_3-hour_Gpi_Cal_3B_42/index.html" + nline
      filAtt_at_Conventions = "None"
      filAtt_at_creation_date = systemfunc("date")
      fileattdef( fnc, filAtt ) ; copy file attributes
    ;===================================================================
    ; predefine the coordinate variables and their dimensionality
    ; note: to get an UNLIMITED record dimension, we set the dimensionality
    ; to -1 and set the unlimited array to True.
    ;===================================================================
        dimNames = (/"time", "lat", "lon"/)
        dimSizes = (/ -1 , nlat, mlon/)
        dimUnlim = (/ True , False, False/)
        filedimdef(fnc,dimNames,dimSizes,dimUnlim)

    ;===================================================================
    ; predefine the the dimensionality of the variables to be written out
    ; filevardef(fnc, "time" ,typeof(time),getvardims(time))
    ;===================================================================
       filevardef(fnc, "time" ,typeof(time), getvardims(time))
       filevardef(fnc, "lat" ,typeof(lat) , getvardims(lat ))
       filevardef(fnc, "lon" ,typeof(lon) , getvardims(lon ))
       filevardef(fnc, "date" ,typeof(date), getvardims(date))
       filevardef(fnc, "datesec" ,typeof(datesec), getvardims(datesec))
       filevardef(fnc, "precip" ,typeof(p) , getvardims(p ))

    ;===================================================================
    ; Copy attributes associated with each variable to the file
    ; All attributes associated with each variable will be copied.
    ;====================================================================
       filevarattdef(fnc,"time",time) ; copy time attributes
       filevarattdef(fnc,"lat",lat) ; copy lat attributes
       filevarattdef(fnc,"lon",lon) ; copy lon attributes
       filevarattdef(fnc,"date",date)
       filevarattdef(fnc,"datesec",datesec)
       filevarattdef(fnc,"precip",p) ; copy p attributes

    ;===================================================================
    ; explicitly exit file definition mode. **NOT REQUIRED**
    ;===================================================================
        setfileoption(fnc,"DefineMode",False)

    ;===================================================================
    ; output only the data values since the dimensionality and such have
    ; been predefined. The "(/", "/)" syntax tells NCL to only output the
    ; data values to the predefined locations on the file.
    ;====================================================================
       fnc->time = (/time/)
       fnc->lat = (/lat /)
       fnc->lon = (/lon /)
       fnc->date = (/date/)
       fnc->datesec= (/datesec/)
       fnc->precip = (/p/)

       delete(P)
       delete(p)
       delete(time)
   end do

end ; only if begin present
Received on Fri Jan 12 2007 - 15:12:03 MST

This archive was generated by hypermail 2.2.0 : Tue Jan 16 2007 - 14:49:45 MST