fils = systemfunc ("ls /backup/seti/setii/NCL-files/1999-2015_DUST04.2003*.nc") ; file paths print(fils) a = addfiles (fils, "r") printVarSummary (a) ListSetType (a, "cat") ; concatenate (=default) ;================================================================================= fils = systemfunc ("ls /backup/seti/setii/NCL-files/1999-2015_SRF.2003*.nc ") ; file paths print(fils) b = addfiles (fils, "r") printVarSummary (b) ListSetType (b, "cat") ; concatenate (=default) ;================================================================================ ; yyyymmddhh ymdhStrt = 2003010100 ; climatology start year/month ;;ymdhLast = 2014123118 ; " last year/month a time = f- ymdhLast = 2003123100 ; " last year/month a time = f- timebnd = a[:]->time_bnds time0 = timebnd(:,0) ; extract the start bound yyyymmddhh= cd_calendar(time0, -3) ;;print(yyyymmddhh) ntStrt = ind(yyyymmddhh.eq.ymdhStrt) ; start time index ntLast = ind(yyyymmddhh.eq.ymdhLast) ; last time index if (ismissing(ntStrt) .or. ismissing(ntStrt)) then print("One or both of the following is missing: ntStrt="+ntStrt+" ntLast="+ntLast) print("Look at yout Start and/or ymdhLast times") exit end if x = a[:]->emflx ; read emission flux from all files x&time = time0 ; override original 'tag' time printVarSummary (x) printMinMax (x,0) print("--------") dimx = dimsizes(x) ntim = dimx(0) nlat = dimx(1) mlon = dimx(2) ;================================================================================ u = b[:]->uas(:,0,:,:) ; read U from all files u&time = time0 ; override original 'tag' time v = b[:]->vas(:,0,:,:) ; read V from all files v&time = time0 ; override original 'tag' time printVarSummary (u) printMinMax (u,0) print("--------") printVarSummary (v) printMinMax (v,0) print("--------") ;===================================================================================== ;;opt = True ;;opt@nval_crit = 4 ; require at least 4 values per day for a daily mean to be calculated. ;;xDay = calculate_daily_values (x, "avg", 0, opt) ;;printVarSummary(xDay) ;;printMinMax (x,0) ;;print("--------") ;; ;;uDay = calculate_daily_values (u, "avg", 0, opt) ;;printVarSummary(uDay) ;;printMinMax (uDay,0) ;;print("--------") ;; ;;vDay = calculate_daily_values (v, "avg", 0, opt) ;;printVarSummary(vDay) ;;printMinMax (vDay,0) ;;print("--------") ;====================================================================================== ; Compute monthly average using all daily 6-hrly data ;====================================================================================== xMonth = calculate_monthly_values(x, "avg", 0, False) printVarSummary (xMonth) printMinMax (xMonth,0) print("--------") uMonth = calculate_monthly_values(u, "avg", 0, False) printVarSummary (uMonth) printMinMax (uMonth,0) print("--------") vMonth = calculate_monthly_values(v, "avg", 0, False) printVarSummary (vMonth) printMinMax (vMonth,0) print("--------") ;************************************************ ; Compute monthly climatology from all year-month means ;************************************************ xClm = clmMonTLL(xMonth) printVarSummary (xClm) ; (12,lat,lon) printMinMax (xClm,0) print("--------") uClm = clmMonTLL(uMonth) printVarSummary (uClm) ; (12,lat,lon) printMinMax (uClm,0) print("--------") vClm = clmMonTLL(vMonth) printVarSummary (vClm) ; (12,lat,lon) printMinMax (vClm,0) print("--------") ;************************************************ ; Explore [statistical dispersion] the monthly climatology ; Also, this information can facilitate setting contour plot limits ;************************************************ opt_sdClm = True opt_sdClm@PrintStat = True stat_x = stat_dispersion(xClm, opt_sdClm ) print("--------") stat_u = stat_dispersion(uClm, opt_sdClm ) print("--------") stat_v = stat_dispersion(vClm, opt_sdClm ) print("--------") ;************************************************ ; Plot ;************************************************* ; NCL graphics will use attached coordinate variables to plot. ; This is not the desired behavior for this dataset. Hence, delete them. ; The xlat, xlon are the lat/lon values for plottong ;************************************************* xlat = a[0]->xlat xlon = a[0]->xlon delete([/xClm&jx,xClm&iy, uClm&jx,uClm&iy, vClm&jx,vClm&iy \ ,xlat&jx, xlat&iy, xlon&jx, xlon&iy /]) latS = min(xlat) latN = max(xlat) lonL = min(xlon) lonR = max(xlon) print("latN="+latN+" latS="+latS) print("lonL="+lonL+" lonR="+lonR) ;************************************************* months = (/"January", "February", "March", "April" \ ,"May", "June", "July", "August" \ ,"September", "October", "November" \ ,"December" /) nmos = dimsizes(months) yyyyStrt = ymdhStrt/1000000 ; climatology start year yyyyLast = ymdhLast/1000000 ; climatology start year ;************************************************* wks = gsn_open_wks("png","monthly_climatology."+yyyyStrt+"-"+yyyyLast) ; send graphics to PNG file ;************************************************* PANEL = True ; draw as panel plot = new (nmos, graphic) ; create graphical array resP = True ; panel options resP@gsnMaximize = True ; maximize image resP@gsnPanelLabelBar = True ; add common colorbar resP@gsnPanelMainString = xClm@long_name+": ["+xClm@units+"]: "+yyyyStrt+"-"+yyyyLast res = True ; plot options desired if (PANEL) then res@gsnDraw = False res@gsnFrame = False else res@tiMainString = xClm@long_name+": ["+xClm@units+"]: "+yyyyStrt+"-"+yyyyLast end if res@gsnAddCyclic = False ; data are not cyclic res@cnFillOn = True ; turn on color fill res@cnInfoLabelOn = False ; turn off contour info label res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off line labels res@cnFillMode = "RasterFill" ; Raster Mode res@cnFillPalette = "WhiteBlueGreenYellowRed" res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = (/5,25,50,75,100,200,300,400,500,750 \ ,1000,1500,2000,2500,3000,3500/) res@lbLabelBarOn = False ; No single label bar res@sfXArray = xlon ; equivalently x@lon2d = xlon res@sfYArray = xlat ; equivalently x@lat2d = xlat res@mpFillOn = False ; turn off gray continents res@mpProjection = "LambertConformal" res@mpLambertParallel1F = a[0]@standard_parallel(0) ; get from the first file [0] res@mpLambertParallel2F = a[0]@standard_parallel(1) res@mpLambertMeridianF = a[0]@longitude_of_projection_origin res@mpLimitMode = "Corners" ; choose region of map res@mpLeftCornerLatF = xlat(0,0) res@mpLeftCornerLonF = xlon(0,0) res@mpRightCornerLatF = xlat(nlat-1,mlon-1) res@mpRightCornerLonF = xlon(nlat-1,mlon-1) res@mpDataBaseVersion = "MediumRes" ; better map outlines res@pmTickMarkDisplayMode= "Always" ; turn on tickmarks ;res@mpOutlineOn = True ;res@mpOutlineBoundarySets = "AllBoundaries" ;res@mpDataBaseVersion = "MediumRes" ; necessary for mpDataSetName to be effective res@gsnLeftString = "" res@gsnRightString = "" resv = True resv@gsnDraw = False ; turn off draw and frame resv@gsnFrame = False ; b/c this is an overlay plot resv@vcRefMagnitudeF = 5.0 ; define vector ref mag resv@vcRefLengthF = 0.045 ; define length of vec ref resv@vcRefAnnoOrthogonalPosF = -0.5 ; move ref vector resv@vcRefAnnoArrowLineColor = "black" ; change ref vector color resv@vcRefAnnoArrowUseVecColor = False ; don't use vec color for ref resv@vcGlyphStyle = "LineArrow" ; turn on curly vectors ;resv@vcLineArrowColor = "white" ; change vector color resv@vcLineArrowColor = "black" ; change vector color ;resv@vcLineArrowThicknessF = 2.0 ; change vector thickness ;resv@vcVectorDrawOrder = "PostDraw" resv@vfXArray = xlon resv@vfYArray = xlat resv = True resv@gsnDraw = False ; turn off draw and frame resv@gsnFrame = False ; b/c this is an overlay plot resv@vcRefMagnitudeF = 5.0 ; define vector ref mag resv@vcRefLengthF = 0.045 ; define length of vec ref resv@vcRefAnnoOrthogonalPosF = -0.5 ; move ref vector resv@vcGlyphStyle = "LineArrow" ; turn on curly vectors ;resv@vcLineArrowColor = "white" ; change vector color resv@vcLineArrowColor = "black" ; change vector color ;resv@vcLineArrowThicknessF = 2.0 ; change vector thickness ;resv@vcVectorDrawOrder = "PostDraw" resv@vfXArray = xlon resv@vfYArray = xlat do nmo=0,nmos-1 ; loop over the months res@gsnCenterString = months(nmo) plot1(nmo) = gsn_csm_contour_map(wks,xClm(nmo,:,:),res) ; create plot plot2(nmo) = gsn_csm_vector(wks,uClm(nmo,:,:),vClm(nmo,:,:),resv) overlay(plot1(nmo),plot2(nmo)) end do shapefile_name1 = "gadm36_IRN_0.shp" shapefile_name2 = "gadm36_IRN_1.shp" lnres = True lnres@gsLineThicknessF = 1.0 lnres@gsLineColor = "Black" shape1 = gsn_add_shapefile_polylines(wks,plot,shapefile_name1,lnres) shape2 = gsn_add_shapefile_polylines(wks,plot,shapefile_name2,lnres) gsn_panel(wks,plot1,(/3,4/),resP) ; now draw as one plot ; create