: This does one day at a time. ;******************************************************************** ; Utility functions to get (1) date information (2) facilitate netCDF ;******************************************************************** function parseFileName ( fNam:string ) ; parse the file names and extract data information ; 1 2a ; 0123456789012345678901234 ; 2A25.20010101.17823.7.HDF local onNam, n, output, tmp_c begin nNam = dimsizes( fNam ) output = new( 3,integer,"No_FillValue") tmp_c = stringtochar(fNam) ; character to integer output(0) = stringtointeger((/tmp_c(5:8)/)) ; YYYY output(1) = stringtointeger((/tmp_c(9:10)/)) ; MM output(2) = stringtointeger((/tmp_c(11:12)/)) ; DD return (output) end ;--------------------------------------------------- function create3d (x[*][*]:numeric, time[*]:numeric) ; convert a 2-dimensional to a 3-dimensional variable ; with time as a coordinate variable. local dimx, ntim, nlat, mlon, x3d begin dimx = dimsizes(x) ntim = dimsizes(time) nlat = dimx(0) mlon = dimx(1) ; perform array expansion x3d = conform_dims( (/ntim,nlat,mlon/), x, (/1,2/)) copy_VarAtts(x, x3d) x3d!0 = "time" x3d&time= time copy_VarCoords(x, x3d(0,:,:) ) ;;print("== CREATE3D ======") ;;printVarSummary(x3d) ; test ;;printMinMax(x3d,0) print("==================") return(x3d) end begin ;************************** ; Main ;************************** PLOT = True netCDF= True nlat = 90 ; 360 mlon = 180 ; 720 lat = latGlobeFo(nlat,"lat","latitude","degrees_north") lon = lonGlobeFo(mlon,"lon","longitude","degrees_east") lon = (/lon -180./) lon&lon=lon ;;print(lat) ;;print(lon) ;***************************************************************** ; Variables to hold binned quantities ;***************************************************************** gbin = new ( (/nlat,mlon/), float ) gknt = new ( (/nlat,mlon/), integer) gbin = 0.0 ; initialization gknt = 0 ;***************************************************************** diri = "/Users/shea/Data/HDF4/" fili = systemfunc("cd "+diri+" ; ls 2A25.20010101*.HDF") ; USER put in day nfil = dimsizes( fili ) ;***************************************************************** ; Loop over the files for the current day only [one day] ;***************************************************************** tStrt = systemfunc("date") ; time the loop (wall clock) vNam = "rain" do nf=0,nfil-1 print(nf+" "+fili(nf)) f = addfile(diri+fili(nf), "r") ; read data xs := f->$vNam$ ; (nscan,nray,ncell1)=>(0,1,2) ; 9239 49 80 x := where(xs.ge.toshort(0), xs*0.01, -9999.0) x@_FillValue = -9999.0 printMinMax(x,0) lat2d := ( f->Latitude ) ; (nscan,nray) lon2d := ( f->Longitude) ; (nscan,nray) LAT1D := ndtooned( conform(x,lat2d, (/0,1/)) ) LON1D := ndtooned( conform(x,lon2d, (/0,1/)) ) print(nf+": "+fili(nf)+": min(lat2d)="+min(lat2d)+" max(lat2d)="+max(lat2d) \ +": min(lon2d)="+min(lon2d)+" max(lon2d)="+max(lon2d)) bin_sum(gbin,gknt,lon,lat,LON1D,LAT1D,ndtooned(x)) end do wallClockElapseTime(tStrt, "Main File Loop", 0) ;***************************************************************** ; User must perform averaging ;***************************************************************** gknt = where(gknt.eq.0 , gknt@_FillValue, gknt) ; avoid division by zero gbin = gbin/gknt ; original units gbin!0 = "lat" gbin!1 = "lon" gbin&lat = lat gbin&lon = lon copy_VarCoords(gbin, gknt) ; copy coords if (isfilevaratt(f, vNam, "long_name")) then gbin@long_name = "BINNED: "+vNam gknt@long_name = "BINNED COUNT: "+vNam end if if (isfilevaratt(f, vNam, "units")) then gbin@units = f->$vNam$@units end if ; change units from mm/hr to mm/day ; (mm/hr)*(24hr/day) = mm/day gbin = gbin*24 gbin@units = "mm/day" ;file_out="out.bin" ; binary output file ;printvarSummary(gbin) ;fbinwrite(file_out,gbin) ;******************************** ;Plot ;******************************** if (PLOT) then GLOBE = False plot = new ( 2, "graphic") pltType = "png" ;pltName = "2A25."+date pltName = "2A25.TEST" wks = gsn_open_wks(pltType ,pltName) gsn_define_colormap(wks, "BlAqGrYeOrReVi200") res = True ; Use plot options res@gsnDraw = False res@gsnFrame = False res@cnFillOn = True ; Turn on color fill res@cnFillMode = "CellFill" ; "RasterFill" ; Turn on raster color res@cnLinesOn = False ; Turn off contourn lines res@mpCenterLonF = 180 ; set map center at 180 res@lbOrientation = "vertical" if (GLOBE) then ; Globe print("********") printMinMax(gbin,0) printMinMax(gknt,0) print("********") res@gsnAddCyclic = True plot(0) = gsn_csm_contour_map_ce(wks,gbin, res) ; globe plot(1) = gsn_csm_contour_map_ce(wks,gknt, res) else ; Region latS = 5 latN = 20 lonL = 115 lonR = 150 res@mpMinLonF = lonL res@mpMaxLonF = lonR res@mpMinLatF = latS res@mpMaxLatF = latN res@gsnAddCyclic = False plot(0) = gsn_csm_contour_map_ce(wks,gbin({latS:latN},{lonL:lonR}), res) plot(1) = gsn_csm_contour_map_ce(wks,gknt({latS:latN},{lonL:lonR}), res) end if ;***************************************************************** ; Panel ;***************************************************************** resP = True ;resP@txString = "TRMM: 2A25: "+yyyy+"/"+month+"/"+day res@gsnMaximize = True ; max size resP@gsnPaperOrientation = "portrait" gsn_panel(wks,plot,(/2,1/), resP) end if exit ;***************************************************************** ; time/date ;***************************************************************** ydhm = parseFileName(fili(0)) yyyy = ydhm(0) month = ydhm(1) day = ydhm(2) hh = 0 mn = 0 ;***************************************************************** ; netCDF ;***************************************************************** if (netCDF) then sec = 0.0 ; COARDS/CF time tunits = "hours since 1990-1-1 00:00:0.0" time = cd_inv_calendar(yyyy,month,day,hh,mn,sec,tunits, 0) time!0 = "time" time@units = tunits time&time = time ;;print(time) ; gregorian date date = yyyy + month + day date@units = "yyyymmdd" date!0 = "time" date&time = time ;;print(date) diro = "./" ; output directory filo = "TRMM2A25Bin_"+yyyy+day+"_bin_sum" ; output file name fout = diro+filo+".nc" system("/bin/rm -f "+fout) ; remove any pre-existing file ncdf = addfile(fout ,"c") ; open output netCDF file ;================================================================ ; create global attributes of the file (not required) ;================================================================ fAtt = True ; assign file attributes fAtt@title = "2A25 Swath Binned: " + yyyy+day fAtt@source_file = "TRMM2A25."+yyyy+day+".nc" ;fAtt@Conventions = "CF 1.0" fAtt@creation_date = systemfunc ("date") fileattdef( ncdf, fAtt ) ; copy file attributes ; recommended (generally) filedimdef(ncdf,"time",-1,True) ; make time UNLIMITED dimension ncdf->time = time ncdf->date = date ncdf->$vNam$ = create3d( gbin, time) kNam = vNam+"_knt" ncdf->$kNam$ = create3d( gknt, time) end if end