;*************************************************************** ; 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 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 ;***************************************************************** ; binned grid definition ; The grid must be constant spacing in lat and in lon. ; The spacing in the lat and lon directions may be different. ; so: 2x2, 3x5, etc ;***************************************************************** nlat = 10800 ; grid coordinates mlon = 10800 lat = latGlobeFo(nlat,"lat","latitude","degrees_north") lon = lonGlobeFo(mlon,"lon","longitude","degrees_east") lon = (/ lon - 180. /) ; subtract 180 from all values lon&lon = lon ; update coordinates ;***************************************************************** ; 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 ;***************************************************************** diri = "./" ; fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2019126.h12v04.006.2019128034342.hdf") fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2002*.hdf") nfil = dimsizes( fili ) print("nfil="+nfil) f0 = addfile(diri+fili(0)+".he2", "r") ; .he2 causes NCL to add lat/lon arrays x = short2flt( f0->$vNam$ ) ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km ) ;***************************************************************** ; Loop over the files: Sum the quantities ; *NO* data filtering is performed: No range checking nor use of AOD_QA_grid1km ;***************************************************************** tStrt = systemfunc("date") ; time the loop (wall clock) do nf=0,nfil-1 ; loop over all files print(nf+" "+fili(nf)) f = addfile(diri+fili(nf)+".he2", "r") ; .he2 causes NCL to add lat/lon arrays ; read data lat2d = f->GridLat_grid1km lon2d = f->GridLon_grid1km x = short2flt( f->$vNam$ ) ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km ) printMinMax(x,0) 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 delete( x ) ; size may change end do wallClockElapseTime(tStrt, "Main Sum Loop: nlat="+nlat+", mlon="+mlon, 0) ;***************************************************************** ; Perform averaging ;***************************************************************** ; 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 ; ;***************************************************************** ;Attempting to print out a netcdf file of the averaged MODIS AOD data for a given time frame ; varName = "Optical_Depth_047_grid1km" varName = GBIN@units ; latloc = 40.821 ; can be more than on location ; lonloc = -73.949 N = 90000000 ; some big number xtxt = new( N,"string","missing") nl = -1 ; do nf=0,nhdf-1 ; f = addfile(pth_hdf(nf)+".he2","r") ; .he2 read HDFEOS information lat2d = f->GridLat_grid1km lon2d = f->GridLon_grid1km ; Using the min max to cover only New York State ; minlat = min(lat2d) ; maxlat = max(lat2d) ; minlon = min(lon2d) ; maxlon = max(lon2d) minlat = 40 maxlat = 46 minlon = -80. maxlon = -71. ; is desired location if (latloc.ge.minlat .and. latloc.le.maxlat .and. \ lonloc.ge.minlon .and. lonloc.le.maxlon ) then x := short2flt(f->$varName$) ; use := syntax in case size changes dimx = dimsizes(x) yyyy = str_get_cols(fil_hdf(nf), 9,12) ddd = str_get_cols(fil_hdf(nf),13,15) norb = dimx(0) xval := rcm2points (lat2d, lon2d, x, latloc, lonloc, 0) nl = nl+1 xtxt(nl) = yyyy+" "+ddd+" "+sprinti("%0.2i", no)+" "+sprintf("%6.2f", xval(no,0)) ; end do ; no if (nl.eq.(N-1)) then print("Array bound of xtxt exceeded: N="+N+" must be made larger") end if end if end do ; nf diro = "./" ; Output directory filo = "test.nc" ; Output file system("/bin/rm -f " + diro + filo) ; remove if exists fout = addfile (diro + filo, "c") ; open output file setfileoption(fout,"DefineMode",True) dimNames = (/"lat", "lon", "GBIN"/) dimSizes = (/ -1 , nlat, nlon, nGBIN /) dimUnlim = (/ True , False, False, False/) filedimdef(fout,dimNames,dimSizes,dimUnlim) filevardef(fout, "GBIN" ,typeof(GBIN),getvardims(GBIN) ) filevardef(fout, "lat" ,typeof(lat),getvardims(lat)) filevardef(fout, "lon" ,typeof(lon),getvardims(lon)) filevarattdef(fout,"GBIN" ,GBIN) ; copy GBIN attributes filevarattdef(fout,"lat" ,lat) ; copy lat attributes filevarattdef(fout,"lon" ,lon) ; copy lon attributes setfileoption(fout,"DefineMode",False) fout->GBIN = (/GBIN/) fout->lat = (/lat/) fout->lon = (/lon/) ;***************************************************************** ; time/date ;***************************************************************** dateinfo= parseFileName(fili(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_binning_full"+yyyymmdd); send graphics to PNG file res = True ; Use plot options res@gsnDraw = False res@gsnFrame = 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 = 180 ; set map center at 180 res@lbOrientation = "vertical" res@mpDataSetName = "Earth..4" ; This database contains ; divisions for other countries. res@mpDataBaseVersion = "MediumRes" ; The following 5 lines are used to create a zoomed image. ; Set limits of map, based on the min/max of the dataset latitude/longitude ; res@mpLimitMode = "LatLon" ; res@mpMinLatF = min(data@lat2d) ; res@mpMaxLatF = max(data@lat2d) ; res@mpMinLonF = min(data@lon2d) ; res@mpMaxLonF = max(data@lon2d) res@mpLimitMode = "LatLon" res@mpMinLatF = 40 res@mpMaxLatF = 46 res@mpMinLonF = -80. res@mpMaxLonF = -71. ; 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.60 res@cnLevelSpacingF = .05 plot(0) = gsn_csm_contour_map(wks,GBIN, res) ;......................................................................... res@cnLevelSelectionMode = "AutoLevels" ; res@cnMinLevelValF = 0. ; res@cnMaxLevelValF = 1. ; res@cnLevelSpacingF = .1 ; plot(1) = gsn_csm_contour_map(wks,GKNT, res) ;.............................................................................. ; plot(0) = gsn_csm_contour_map(wks,GBIN, res) plot(1) = gsn_csm_contour_map(wks,GKNT, res) ;***************************************************************** ; Panel ;***************************************************************** dlat = lat(1)-lat(0) resP = True resP@gsnPanelMainString = "MCD19A2: "+yyyy+"/"+month+"/"+day+ \ ": "+nlat+"x"+mlon+" dlat="+sprintf("%5.3f", dlat) resP@gsnMaximize = True ; max size ;resP@gsnPaperOrientation = "portrait" gsn_panel(wks,plot,(/2,1/), resP) end if