; Ozone site plotting for R values of DEC verses GEOS-Chem scatter plot values to display the spatial difference across NYS. ; ; ; ;---------------------------------------------------------------------- ; Procedure for adding a labelbar at a given NDC location, given ; the levels and colors to use. ;---------------------------------------------------------------------- undef("add_labelbar") procedure add_labelbar(wks,plot,levels,colors) local lbres, labels begin nlevels = dimsizes(levels) ;---------------------------------------------------------------------- ; Draw a labelbar ;---------------------------------------------------------------------- lbres = True lbres@vpWidthF = 0.80 ; width lbres@vpHeightF = 0.10 ; height lbres@lbPerimOn = False ; Turn off perimeter. lbres@lbOrientation = "Horizontal" ; Default is vertical. lbres@lbLabelAlignment = "InteriorEdges" ; Default is "BoxCenters" lbres@lbFillColors = colors ; Colors for boxes. lbres@lbMonoFillPattern = True ; Fill them all solid. lbres@lbLabelFontHeightF = 0.012 ; label font height labels = sprintf("%4.2f",levels) lbid = gsn_create_labelbar(wks,nlevels+1,labels,lbres) ; ; Now, create some annotation resources indicating how we want to ; attach the labelbar to the plot. Here, we are using the top right ; corner of the labelbar as the point which we are going to position ; it, and then we use amParallelPosF and amOrthogonalPosF to indicate ; where we want to place it. ; ; amParallelPosF/amOrthogonalPosF ; ; 0.0/ 0.0 - annotation in dead center of plot ; 0.5/ 0.5 - annotation at bottom right of plot ; 0.5/-0.5 - annotation at top right of plot ; -0.5/-0.5 - annotation at top left of plot ; -0.5/ 0.5 - annotation at bottom left of plot ; amres = True amres@amJust = "TopCenter" amres@amParallelPosF = 0.0 ; keep labelbar centered amres@amOrthogonalPosF = 0.6 ; move down and outside of plot ; ; Give both annotation id and labelbar id unique names. ; ; Attaching them to plot with unique names ensures that ; labelbar "lives" outside this procedure. ; tmpid1 = "anno"+unique_string("id") tmpid2 = "lbar"+unique_string("id") plot@$tmpid1$ = gsn_add_annotation(plot,lbid,amres) plot@$tmpid2$ = lbid end ;---------------------------------------------------------------------- ; Main code. ;---------------------------------------------------------------------- begin ;--------------------------- ; Loading input data file strs = asciiread("MFB-O3Montlhy-APR.txt",-1,"string") delim = " " nfields = str_fields_count(strs(0), delim) ; print(nfields) fieldavg = 3 ;avged data R = tofloat(str_get_field(strs,3," ")) R@_FillValue=integertoshort(-999) print(R) lat = tofloat(str_get_field(strs,1," ")) print(lat) Nlat = lat lon = tofloat(str_get_field(strs,2," ")) print(lon) Nlon = lon ;--------------------------- ; Set up markers ; levels = (/0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/) levels = (/-40.0,-20.0,-10.0,-7.0,-5.0,-2.0,-1.0,0.0,1.0,4.0/) nlevels = dimsizes(levels) ; colors = read_colormap_file("WhViBlGrYeOrRe") ; setting color markers colors = (/"10","30","38","48","56","66","74","94","450","600"/) labels = new(dimsizes(levels)+1,string) ; Labels for legend. ;------------------------------------------------------------- num_distinct_markers = dimsizes(levels)+1 ; number of distinct markers lat_new = new((/num_distinct_markers,dimsizes(R)/),float,-999) lon_new = new((/num_distinct_markers,dimsizes(R)/),float,-999) ; print(arr) ; print(lat_new) ; print(lon_new) ; ; Group the points according to which range they fall in. At the ; same time, create the label that we will use later in the legend. ; do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(R.lt.levels(i)) labels(i) = "x < " + levels(i) end if if (i.eq.num_distinct_markers-1) then indexes = ind(R.ge.max(levels)) labels(i) = "x >= " + max(levels) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(R.ge.levels(i-1).and.R.lt.levels(i)) labels(i) = levels(i-1) + " <= x < " + levels(i) end if ; ; Now that we have the set of indexes whose values fall within ; the given range, take the corresponding lat/lon values and store ; them, so later we can color this set of markers with the appropriate ; color. ; if (.not.any(ismissing(indexes))) then npts_range = dimsizes(indexes) ; # of points in this range. lat_new(i,0:npts_range-1) = lat(indexes) lon_new(i,0:npts_range-1) = lon(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do ;============================================================ ; Start the graphics. wks = gsn_open_wks("png","MFB-O3Montlhy-APR.png") ; send graphics to PNG file ; gsn_define_colormap(wks,"gui_default") ; define a different color map gsn_define_colormap(wks,"WhViBlGrYeOrRe") ; ; plot = new(1,graphic) ; Set up some map resources. ; mpres = True mpres@gsnDraw = False ; Will draw later mpres@gsnMaximize = True ; Maximize plot in frame. mpres@gsnFrame = False ; Don't advance the frame mpres@pmTickMarkDisplayMode = "Always" ; We want county outlines. mpres@mpDataBaseVersion = "MediumRes" mpres@mpDataSetName = "Earth..2" mpres@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres@mpOutlineSpecifiers = (/"Land","New York:counties"/) ; Zoom in on NY, U.S. mpres@mpLimitMode = "LatLon" mpres@mpMinLatF = 40 mpres@mpMaxLatF = 46 mpres@mpMinLonF = -80.0 mpres@mpMaxLonF = -71.0 ; Turn on tickmarks. mpres@gsnTickMarksOn = True mpres@mpFillColors = (/"transparent","transparent", \ "lightgray","transparent"/) ;assign light gray to land masses ; Main title mpres@tiMainString = "DEC Ozone Site Winter R Values" map = gsn_csm_map(wks,mpres) ; ; Create logical variables to hold the marker and text resources. ; These markers are different than the XY markers, because they are not ; associated with an XY plot. You can put these markers on any plot. ; mpres = True mpres@gsMarkerIndex = 16 ; Use filled dots for markers. gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. ; txres = True ; txres@txFontHeightF = 0.015 ; txres = True ; txres@txFontHeightF = 0.015 ; ; Loop through each grouping of markers, and draw them one set at ; a time, assigning the proper color and size with gsn_marker. ; ; At the same time, draw a legend showing the meaning of the markers. ; ; xleg = (/0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.8,0.9,0.9/) ; Location of ; xtxt = (/0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.7,0.9,0.9,1.0,1.0/) ; legend markers ; xleg = (/0.1,0.1,0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.8,0.9,0.9/) ; Location of ; xtxt = (/0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5,0.6,0.6,0.7,0.7,0.8,0.7,0.9,0.9,1.0,1.0/) ; 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. ;---------------------------------------------- base_size = 1 pmid = new(num_distinct_markers,graphic) do i = 0, num_distinct_markers-1 if (.not.ismissing(lat_new(i,0))) gsres@gsMarkerColor = colors(i) gsres@gsMarkerThicknessF = 0.7*(i+1) ; gsres@gsMarkerColor = colors(i) ; gsres@gsMarkerThicknessF = 0.7*(i+1) ; gsres@gsMarkerThicknessF = 0.7 pmid(i) = gsn_add_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 labelbar and advance frame. add_labelbar(wks,map,levels,colors) draw(map) frame(wks) ; Advance the frame. end