load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" ;******************************************** begin a = addfile("./model/20260728/gefs.chem.t00z.a2d_0p25.f018.grib2","r") z = addfile("./Airnowobs/20260728/Airnow_HrlyPM25_2026072818.nc","r") ;********************* ;---------------------------------------------------------------------- ; This function attaches a labelbar to the given plot. ;---------------------------------------------------------------------- function attach_labelbar(wks,map,arr[*]:numeric,colors[*]) local lbres, vph, vpw, nboxes begin getvalues map "vpHeightF" : vph "vpWidthF" : vpw end getvalues nboxes = dimsizes(colors) lbres = True ; labelbar only resources lbres@lbAutoManage = False ; Necessary to control sizes lbres@lbFillColors = colors lbres@vpWidthF = 0.65 * vpw ; labelbar width lbres@vpHeightF = 0.2 * vph ; labelbar height lbres@lbMonoFillPattern = True ; Solid fill pattern lbres@lbLabelFontHeightF = 0.013 ; font height. default is small lbres@lbOrientation = "horizontal" lbres@lbPerimOn = False lbres@lbLabelAlignment = "InteriorEdges" lbid = gsn_create_labelbar(wks,nboxes,""+tofloat(arr),lbres) amres = True amres@amJust = "TopCenter" amres@amParallelPosF = 0.0 ; Center amres@amOrthogonalPosF = 0.6 ; Bottom annoid = gsn_add_annotation(map,lbid,amres) return(annoid) end gsd1 = a->AOTK_P48_L10_GLL0_A62000 printVarSummary(gsd1) gsd1=gsd1(::-1,:) ; printVarSummary(gsd1) ; gsd1 = lonFlip(gsd1) printVarSummary(gsd1) pm25_data = z->obs_val hdr_lat = z->hdr_lat hdr_lon = z->hdr_lon obs_hid = z->obs_hid obs_vid = z->obs_vid ; print(obs_vid) new_var = ind(obs_vid .eq. 1) ; print(new_var) obs_lat = (/ hdr_lat(obs_hid(:)) /) ; now indexed by nobs obs_lon = (/ hdr_lon(obs_hid(:)) /) ; now indexed by nobs obs_val = z->obs_val data_pm = (/pm25_data(new_var(:))/) obsv_lat = (/ obs_lat(new_var(:)) /) ; now indexed by nobs obsv_lon = (/ obs_lon(new_var(:)) /) printVarSummary(obsv_lat) printVarSummary(obsv_lon) printVarSummary(pm25_data) printVarSummary(data_pm) arr = (/5.,15.,25.,35.,50.,60.,80.,100.,120.,150.,200./) ; bin settings (bin0 = < 0., ; bin1 = 0.:4.999, etc.) colors = (/10,30,38,48,56,66,74,94,114,134,154,174/) ; marker colors, dimsizes must ; be equal to dimsizes(arr)+1 labels = new(dimsizes(arr)+1,string) ; Labels for legend. narr = dimsizes(arr) num_distinct_markers = dimsizes(arr)+1 ; number of distinct markers lat_new = new((/num_distinct_markers,dimsizes(data_pm)/),float,-9999) lon_new = new((/num_distinct_markers,dimsizes(data_pm)/),float,-9999) printVarSummary(num_distinct_markers) do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(data_pm.lt.arr(5)) labels(i) = "x < " + arr(5) end if if (i.eq.num_distinct_markers-1) then indexes = ind(data_pm.ge.max(arr)) labels(i) = "x >= " + max(arr) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(data_pm.ge.arr(i-1).and.data_pm.lt.arr(i)) labels(i) = arr(i-1) + " <= x < " + arr(i) end if if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) ; # of points in this range. printMinMax (indexes, 0) print ("num missing = " + num (ismissing (indexes))) print (i) print (npts_range) printVarSummary(lat_new) ; printVarSummary(obs_lat) printVarSummary(indexes) lat_new(i,0:npts_range-1) = obs_lat(indexes) lon_new(i,0:npts_range-1) = obs_lon(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do wks_type = "png" wks_type@wkWidth = 3500 wks_type@wkHeight = 3500 wname = "Testfig_PM25_20260728" ;; wks = gsn_open_wks("x11","png") wks = gsn_open_wks(wks_type,wname) gsn_define_colormap(wks,"percent_11lev") ;************************************************ ; create panel plots ;************************************************* plot_A = new (1, graphic) ; create graphical array plot_B = new (1, graphic) mpres = True mpres@gsnMaximize = True ; Maximize plot in frame. mpres@gsnFrame = False ; Don't advance the frame mpres@mpFillDrawOrder = "PostDraw" ; Zoom in on United States. mpres@mpMinLonF = 0.0 ; set your minimum longitude mpres@mpMaxLonF = 360.0 ; set your maximum longitude (> 180) mpres@mpCenterLonF = 180.0 ; set the center to match the middle (e.g., 180) mpres@mpMinLatF = 20. mpres@mpMaxLatF = 65. mpres@mpMinLonF = 220. mpres@mpMaxLonF = 300. mpres@mpFillColors = (/"transparent","transparent","gray","transparent"/) mpres@mpFillOn = False map = gsn_csm_map(wks,mpres) ; contour = gsn_csm_contour(wks,gsd1({25:55},{235:290}),cnres) ; overlay(map,contour) ; mpres@gsnFrame = False gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. gsres@gsMarkerSizeF = 0.007 ;---Get nice spacing through color map for marker colors getvalues wks "wkColorMapLen" : clen ; number of colors in color map end getvalues nstep = (clen-2)/narr colors = ispan(2,clen-1,nstep) ;---Loop through each "bin" and attach the markers to the map. do i = 0, num_distinct_markers-1 if (.not.ismissing(lat_new(i,0))) gsres@gsMarkerColor = colors(i) dumstr = unique_string("marker") map@$dumstr$ = gsn_add_polymarker(wks,map(0),lon_new(i,:),lat_new(i,:),gsres) end if end do gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. txres = True txres@txFontHeightF = 0.015 do i = 0, num_distinct_markers-1 if (.not.ismissing(lat_new(i,0))) gsres@gsMarkerColor = colors(i) dumstr = unique_string("marker") map@$dumstr$ = gsn_add_polymarker(wks,map(0),lon_new(i,:),lat_new(i,:),gsres) end if end do lbid = attach_labelbar(wks,map,arr,colors) ; Attach a labelbar ; xleg = (/0.07,0.07,0.32,0.32,0.57,0.57,0.82,0.82,0.93,0.93,0.97,0.97/) ; Location of ; xtxt = (/0.16,0.16,0.41,0.41,0.66,0.66,0.91,0.91,0.95,0.95,0.99,0.99/) ; legend markers ; yleg = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/) ; and text ; ytxt = (/0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17,0.22,0.17/) ; strings. ; do i = 0, num_distinct_markers-1 ; if (.not.ismissing(lat_new(i,0))) ; gsres@gsMarkerColor = colors(i) ; gsres@gsMarkerThicknessF = 0.7*(i+1) ; gsn_polymarker(wks,map,lon_new(i,:),lat_new(i,:),gsres) ; ; Add marker and text for the legend. ; ; gsn_polymarker_ndc(wks, xleg(i),yleg(i),gsres) ; gsn_text_ndc (wks,labels(i),xtxt(i),ytxt(i),txres) ; end if ; end do draw(map) frame(wks) delete(wks) system("convert -trim " + wname + ".png " + wname + "_f018.png") ;******************************************** end