;t*************************************************************** ; binning_2.ncl ; ; Concepts illustrated: ; - Reading multiple fileslwith data from MODIS "MCD19A2" ; - Binning and summing the data from each pass [swath] ; - Averaging binned data ; - Plotting ; - Creating a netCDF (or HDF) file ; ;*************************************************************** ; Utility functions to get (1) date information (2) facilitate netCDF ;******************************************************************** ; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf ; YYYYDDD = Year and Day of Year of acquisition ; hHH = Horizontal tile number (0-35) ; vVV = Vertical tile number (0-17) ; CCC = Collection number ; YYYYDDDHHMMSS = Production Date and Time ; Example: ; MCD19A2.A2019161.h32v11.006.2019163201550.hdf ;******************************************************************** undef("parseFileName") function parseFileName ( fNam[*]:string ) ; parse the file name(s) and extract data information ; https://www.ncl.ucar.edu/Document/Functions/Built-in/str_get_cols.shtml ; 012345678901234567890123456789012345678901234 ; 1 2 3 4 ; MCD19A2.AYYYYDDD.hHHvVV.CCC.YYYYDDDHHMMSS.hdf ; YYYYDDD tile local n, nNam, time_info, YYYY, DDD, yyyyddd, yyyymmdd, mmdd, mm, dd begin nNam = dimsizes( fNam ) YYYY = toint(str_get_cols(fNam, 9,12)) ; YYYY DDD = toint(str_get_cols(fNam,13,15)) ; DDD yyyyddd = YYYY*1000 + DDD yyyymmdd = yyyyddd_to_yyyymmdd(yyyyddd) mmdd = yyyymmdd-(yyyymmdd/10000)*10000 mm = mmdd/100 dd = mmdd-mm*100 time_info = new( (/nNam,6/), integer,"No_FillValue") do n=0,nNam-1 time_info(n,0) = yyyyddd time_info(n,1) = YYYY time_info(n,2) = DDD time_info(n,3) = yyyymmdd time_info(n,4) = mm time_info(n,5) = dd end do return (time_info) end ;************************************************************** ; MAIN ;************************************************************** vNam = "Optical_Depth_047_grid1km" PLOT = True latS = 40 latN = 46 lonW = -80 lonE = -71 nlat = toint((latN-latS)/0.01) ; 1km <~> 0.01 deg (approx) mlon = toint((lonE-lonW)/0.01) nlat = 100 ; larger bins mlon = 100 ; nlat and mlon can be different lat = fspan(latS,latN,nlat) lon = fspan(lonW,lonE,mlon) lat!0 = "lat" ; meta data lon!0 = "lon" lat@long_name = "latitude" lon@long_name = "longitude" lat@units = "degrees_north" lon@units = "degrees_east" lat&lat = lat lon&lon = lon printVarSummary(lat) printMinMax(lat,0) print("=====") printVarSummary(lon) printMinMax(lon,0) print("=====") ;***************************************************************** ; Variables to hold binned quantities ;***************************************************************** GBIN = new ( (/nlat,mlon/), float ) GKNT = new ( (/nlat,mlon/), integer ) GBIN = 0.0 ; initialization GKNT = 0 ;***************************************************************** ; File info ; HDFEOSVersion : HDFEOS_V2.19 ; StructMetadata_0 : GROUP=SwathStructure ;***************************************************************** ; Examine the 1st file for assorted information ;***************************************************************** dirMCD = "/Users/shea/Data/HDF4/" ; filMCD = systemfunc("cd "+dirMCD+" ; ls MCD19A2.A2019126.h12v04.006.2019128034342.hdf") filMCD = systemfunc("cd "+dirMCD+" ; ls MCD19A2.A2002*.hdf") nfil = dimsizes( filMCD ) print(filMCD) print("nfil="+nfil) print("=========") f0 = addfile(dirMCD+filMCD(0)+".he2", "r") ; .he2 causes NCL to add lat/lon arrays x0 = short2flt( f0->$vNam$ ) ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km ) lat2d0 = f0->GridLat_grid1km lon2d0 = f0->GridLon_grid1km dimx0 = dimsizes(x0) printVarSummary(x0) printMinMax(x0,0) print("=========") printVarSummary(lat2d0) printMinMax(lat2d0,0) print("=========") printVarSummary(lon2d0) printMinMax(lon2d0,0) print("=========") print(dimx0) print("=========") norbit0 = dimx0(0) nYdim = dimx0(1) nXdim = dimx0(2) delete([/x0, lat2d0, lon2d0, dimx0 /]) ; no longer needed ;***************************************************************** ; [1] Loop over all the files ; [2] For each file: ; (a) loop over all the orbits; NOTE: each orbit has a different lat2d/lon2d values ; (b) Sum/Bin the quantities using NCL functions ; [3] *NO* data filtering is performed: No range checking ;***************************************************************** tStrt = systemfunc("date") ; time the loop (wall clock) do nf=0,nfil-1 ; loop over all files print(nf+" "+filMCD(nf)) f = addfile(dirMCD+filMCD(nf)+".he2", "r") ; .he2 causes NCL to add lat/lon arrays ; read data lat2d = f->GridLat_grid1km ; 1200 lon2d = f->GridLon_grid1km ; 1200 x := short2flt( f->$vNam$ ) ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km ) dimx = dimsizes(x) norbit = dimx(0) do no=0,norbit-1 ; loop over all orbits within the current file bin_sum(GBIN,GKNT,lon,lat \ ,ndtooned(lon2d), ndtooned(lat2d),ndtooned(x(no,:,:)) ) end do end do ; nf wallClockElapseTime(tStrt, "Main Sum Loop: nlat="+nlat+", mlon="+mlon, 0) ;***************************************************************** ; Perform averaging: This is a global grid ;***************************************************************** ; avoid division by 0.0 GKNT = where(GKNT.eq.0 , GKNT@_FillValue, GKNT) GBIN = GBIN/GKNT 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 AVERAGE: "+vNam GKNT@long_name = "BINNED COUNT: "+vNam end if if (isfilevaratt(f, vNam, "units")) then GBIN@units = f->$vNam$@units end if ;***************************************************************** ; Examine the global grid ;***************************************************************** printVarSummary(GBIN) printMinMax(GBIN,0) print("--------") opt_sd = True opt_sd@PrintStat = True stat_all = stat_dispersion(GBIN, opt_sd ) print("--------") ;***************************************************************** ; Examine the region of interest ;***************************************************************** printVarSummary(GBIN({latS:latN},{lonW:lonE})) printMinMax(GBIN({latS:latN},{lonW:lonE}),0) stat_regional = stat_dispersion(GBIN({latS:latN},{lonW:lonE}), opt_sd ) print("--------") printVarSummary(GKNT({latS:latN},{lonW:lonE})) printMinMax(GKNT({latS:latN},{lonW:lonE}),0) stat_regional = stat_dispersion(GKNT({latS:latN},{lonW:lonE}), opt_sd ) print("--------") ;***************************************************************** ; Write out a regional grid ;***************************************************************** dirnc = "./" ; Output directory filnc = "test_regional.nc" ; Output file pthnc = dirnc + filnc system("/bin/rm -f " + pthnc) ; remove if exists fnc = addfile (dirnc + filnc, "c") ; open output file fnc->GBIN = GBIN({latS:latN},{lonW:lonE}) dirhdf = "./" ; Output directory filhdf = "test.regional.hdf" ; Output file pthhdf = dirhdf+ filhdf system("/bin/rm -f " + pthhdf) ; remove if exists fhdf = addfile (dirhdf + filhdf, "c") ; open output file fhdf->GBIN = GBIN({latS:latN},{lonW:lonE}) ;***************************************************************** ; time/date ;***************************************************************** dateinfo= parseFileName(filMCD(0)) yyyyddd = dateinfo(0,0) ; this is the date info for the 1st file yyyy = dateinfo(0,1) ddd = dateinfo(0,2) yyyymmdd= dateinfo(0,3) month = dateinfo(0,4) day = dateinfo(0,5) ;***************************************************************** ; PLOTS ;***************************************************************** if (PLOT) then ;; setvalues NhlGetWorkspaceObjectId() ;; "wsMaximumSize" : 900000000 ;; end setvalues plot = new ( 2, "graphic") wks = gsn_open_wks("png", "MCD19A2_binned_region."+yyyymmdd); send graphics to PNG file res = True ; Use plot options res@gsnDraw = False res@gsnFrame = False res@gsnAddCyclic = False res@cnFillOn = True ; Turn on color fill ;;res@cnFillPalette = "wh-bl-gr-ye-re" ; set color map res@cnFillPalette = "BlAqGrYeOrReVi200" ; set color map res@cnFillMode = "RasterFill" ; Turn on raster color res@cnLinesOn = False ; Turn off contourn lines res@mpCenterLonF = 0 res@lbOrientation = "vertical" res@mpDataBaseVersion = "MediumRes" ;;res@mpDataSetName = "Earth..4" ; This database contains ; divisions for other countries. res@mpLimitMode = "LatLon" res@mpMinLatF = latS res@mpMaxLatF = latN res@mpMinLonF = lonW res@mpMaxLonF = lonE ; res@mpUSStateLineColor = "gray" res@mpUSStateLineDashPattern = 7 res@mpUSStateLineThicknessF = 1 ; res@mpNationalLineThicknessF= 2. ; double the thickness of national boundaries res@mpNationalLineThicknessF= 1. ; double the thickness of national boundaries res@mpOutlineBoundarySets = "National" ; draw national boundaries res@mpOutlineSpecifiers = (/"Canada : Provinces","United States : States"/) ; draw US States, Canadian provinces ; res@mpOutlineSpecifiers = (/"China:states","Mongolia","India","Nepal","Bhutan","Bangladesh","Myanmar", \ ; "Thailand","Cambodia","Vietnam","Taiwan","Japan"/) ;......................................................................... res@cnLevelSelectionMode = "ManualLevels" res@cnMinLevelValF = 0. res@cnMaxLevelValF = 0.70 res@cnLevelSpacingF = 0.025 ;;plot(0) = gsn_csm_contour_map(wks,GBIN, res) plot(0) = gsn_csm_contour_map(wks,GBIN({latS:latN},{lonW:lonE}), res) ;......................................................................... ; these are a functo of bin size ;......................................................................... ;;res@cnLevelSelectionMode = "AutoLevels" ;;res@cnMinLevelValF = 4. ;;res@cnMaxLevelValF = 30. ;;res@cnLevelSpacingF= 2. res@cnMinLevelValF = 25. res@cnMaxLevelValF = 130. res@cnLevelSpacingF= 5. ;;plot(1) = gsn_csm_contour_map(wks,GKNT, res) plot(1) = gsn_csm_contour_map(wks,GKNT({latS:latN},{lonW:lonE}), res) ;***************************************************************** ; Panel ;***************************************************************** dlat = lat(1)-lat(0) resP = True resP@gsnMaximize = True ; max size resP@gsnPanelMainString = "MCD19A2: "+yyyy+"/"+month+"/"+day+ \ ": "+nlat+"x"+mlon+" dlat="+sprintf("%5.3f", dlat) ;resP@gsnPaperOrientation = "portrait" gsn_panel(wks,plot,(/2,1/), resP) end if