;*************************************************************** ; 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 ;************************************************************** begin vNam = "Optical_Depth_047_grid1km" PLOT = True ;********** New code addition for printing location data********************** latloc = 40.821 ; can be more than on location lonloc = -73.949 ; N = 20000 ; some big number N = 1 xtxt = new( N,"string","missing") nl = -1 ;**************************************************************************** ;***************************************************************** ; 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 = 1200 ; grid coordinates ; mlon = 1200 nlat = 10800 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 ;***************************************************************** ; creating a new array to place numerous days/months/years together ; X01 = new( (/4,nlat,mlon/),"integer") ; X02 = new((/4,1200,1200/),"integer") GBIN = new ( (/nlat,mlon/), float ) GKNT = new ( (/nlat,mlon/), integer ) ; X01 = new ( (/nlat,mlon/), integer ) GBIN = 0.0 ; initialization GKNT = 0 X01 = 0 ;***************************************************************** ; File info ; HDFEOSVersion : HDFEOS_V2.19 ; StructMetadata_0 : GROUP=SwathStructure ;***************************************************************** diri = "./" ; fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2017259*.hdf") ; fili = systemfunc("cd "+diri+" ; ls MCD19A2.A2017*.hdf") ;************ multiple year averaging **************************** fili = systemfunc("cd "+diri+" ; ls MCD19A2.A*.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 ) printVarSummary(x) printMinMax(x,0) print("=============================") nmsg = num(ismissing(x)) ngood= num(.not.ismissing(x)) npts = nmsg+ngood npc = (ngood/tofloat(npts))*100 print("nmsg ="+nmsg) print("ngood="+ngood) print("npc (% good)="+npc) print("=============================") nx = num(x.gt.-100 .and. x.lt.5000) ; # of values bwtweem -100 and 5000 after scale fac -.01 and 5 ; ***************************************************************** nx = num(x.ge. -0.1 .and. x.le.5.0) ; values after unpacking ; within allowed range print("nx="+nx+": this should match ngood="+ngood) ; ie: there are no out of range values ; X01 = where(x.ge. -0.1 .and. x.le.5.0, 1, 0) printVarSummary(X01) n1n = num(X01.eq.1) n1s = sum(X01) ; X02 = nls print("n1n="+n1n+" n1s="+n1s) print("=============================") ;***************************************************************** ; 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 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( f0->$vNam$ ) ; ( Orbits_grid1km, YDim_grid1km, XDim_grid1km ) printMinMax(x,0) dimx = dimsizes(x) norbit = dimx(0) ;--------------------------------------------------------------------------------------------- GBIN = new ( (/nlat,mlon/), float ) GKNT = new ( (/nlat,mlon/), integer ) ; X01 = new ( (/nlat,mlon/), integer ) GBIN = 0.0 ; initialization GKNT = 0 X01 = 0 ;--------------------------------------------------------------------------------------------- ; do no=0,norbit-1 ; loop over all orbits within the current file ; do no=3,3 ; loop over all orbits within the current file do no=0,norbit-1 bin_sum(GBIN,GKNT,lon,lat \ ,ndtooned(lon2d), ndtooned(lat2d),ndtooned(x(no,:,:)) ) end do delete( x ) ; size may change ; delete( X01 ) ; 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 ;***************************************************************** ;***************************************************************** ;******************** location data grabbing code portition*************** minlat = min(lat2d) maxlat = max(lat2d) minlon = min(lon2d) maxlon = max(lon2d) if (latloc.ge.minlat .and. latloc.le.maxlat .and. \ lonloc.ge.minlon .and. lonloc.le.maxlon ) then x := short2flt(f0->$vNam$) ; use := syntax in case size changes ; x := short2flt(f->$vNam$) ; use := syntax in case size changes dimx = dimsizes(x) yyyy = str_get_cols(fili(nf), 9,12) ddd = str_get_cols(fili(nf),13,15) norb = dimx(0) ; xval := rcm2points (lat2d, lon2d, GBIN, latloc, lonloc, 0) xval := rcm2points (lat2d, lon2d, x, latloc, lonloc, 0) xval := avg(xval) print(xval) ; nl = nl+1 printVarSummary(xval) xtxt = yyyy+" "+ddd+" "+sprinti("%0.2i", no)+" "+sprintf("%6.2f", xval(0)) ; xtxt(nl) = yyyy+" "+ddd+" "+sprinti("%0.2i", no)+" "+sprintf("%6.2f", xval(0)) ; xtxt(nl) = yyyy+" "+ddd+" "+sprinti("%0.2i", no)+" "+sprintf("%6.2f", xval(no)) ; xtxt(nl) = yyyy+" "+ddd+" "+sprinti("%0.2i", no)+" "+sprintf("%6.2f",xval(no,:)) dateinfo= parseFileName(fili(nf)) 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) ;***************** Writing output for site location ****************** dirtxt = "./" ; filtxt = "MCD19A2.OD47.:"+yyyymmdd+ \ ": +round(latloc*10,3)+"_"+abs(round(lonloc*10,3))+".txt" filtxt = "MCD19A2.OD47.Location-CCNY"+yyyymmdd+".txt" pthtxt = dirtxt+filtxt system("/bin/rm -f "+pthtxt) ; remove any preexisting asciiwrite(pthtxt, xtxt(0)) ; asciiwrite(pthtxt, xtxt(0:nl)) end if ;***************************************************************** ;***************************************************************** ; if (isfilevaratt(x, vNam, "long_name")) then if (isfilevaratt(f, vNam, "long_name")) then GBIN@long_name = "BINNED AVERAGE: "+vNam GKNT@long_name = "BINNED COUNT: "+vNam ; X01@long_name = "BINNED COUNT: "+vNam end if if (isfilevaratt(f, vNam, "units")) then GBIN@units = f->$vNam$@units end if printVarSummary(GKNT) ;***************************************************************** ; time/date ;***************************************************************** ; dateinfo= parseFileName(fili(0)) dateinfo= parseFileName(fili(nf)) 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) ;***************** Writing output for site location ****************** ; dirtxt = "./" ; filtxt = "MCD19A2.OD47.:"+yyyymmdd+ \ ": +round(latloc*10,3)+"_"+abs(round(lonloc*10,3))+".txt" ; filtxt = "MCD19A2.OD47.:"+yyyymmdd+".txt" ; pthtxt = dirtxt+filtxt ; system("/bin/rm -f "+pthtxt) ; remove any preexisting ; asciiwrite(pthtxt, xtxt(0)) ; asciiwrite(pthtxt, xtxt(0:nl)) ;***************************************************************** ; PLOTS ;***************************************************************** if (PLOT) then setvalues NhlGetWorkspaceObjectId() "wsMaximumSize" : 1000000000 end setvalues plot = new ( 2, "graphic") wks = gsn_open_wks("png", "MCD19A2_binning_NYS"+yyyymmdd); send graphics to PNG file ; wks = gsn_open_wks("png", "MCD19A2_binning_NYC"+yyyymmdd); send graphics to PNG file ; wks = gsn_open_wks("png", "MCD19A2_binning_ALB"+yyyymmdd); send graphics to PNG file ; wks = gsn_open_wks("png", "MCD19A2_binning_SYR"+yyyymmdd); send graphics to PNG file ; wks = gsn_open_wks("png", "MCD19A2_binning_BGM"+yyyymmdd); send graphics to PNG file ; wks = gsn_open_wks("png", "MCD19A2_binning_BUF"+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@lbLabelAutoStride=True ; ensure labels do not overlap ;......................................................................... ; res@cnLevelSelectionMode = "ManualLevels" ; res@cnMinLevelValF = 0. ; res@cnMaxLevelValF = 0.6 ; 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(0) = gsn_csm_contour_map(wks,GKNT, res) ;.............................................................................. res@cnLinesOn = False ; Turn off contourn lines 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) ;***************** Full NYS Plot *************** res@mpLimitMode = "LatLon" res@mpMinLatF = 40 res@mpMaxLatF = 46 res@mpMinLonF = -80. res@mpMaxLonF = -71. ;************** NYC Plot Only ***************** ; res@mpLimitMode = "LatLon" ; res@mpMinLatF = 40 ; res@mpMaxLatF = 41 ; res@mpMinLonF = -74.5 ; res@mpMaxLonF = -73.5 ; 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@mpCenterLonF = 180 ; set map center at 180 res@lbOrientation = "vertical" ;......................................................................... 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) ; plot(2) = gsn_csm_contour_map(wks, X01, res) ;------------------------------------------------------------------------------------------------------------ delete (GBIN) delete (GKNT) ;------------------------------------------------------------------------------------------------------------ ;***************************************************************** ; 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 end do end