; 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("DEC-GEOS-NOx-FB-Winter.txt",-1,"string") strs1 = asciiread("DEC-GEOS-NOx-FB-Spring.txt",-1,"string") strs2 = asciiread("DEC-GEOS-NOx-FB-Summer.txt",-1,"string") strs3 = asciiread("DEC-GEOS-NOx-FB-Fall.txt",-1,"string") ;Panel 1 Read in data delim = " " nfields = str_fields_count(strs(0), delim) ; print(nfields) fieldavg = 4 ;avged data R = tofloat(str_get_field(strs,4," ")) R@_FillValue=integertoshort(-999) print(R) lat = tofloat(str_get_field(strs,2," ")) print(lat) ; Nlat = new((/51,1/),"float",-999.) Nlat = lat ; print(Nlat) ; lat@_FillValue=integertoshort(-999) lon = tofloat(str_get_field(strs,3," ")) print(lon) ; Nlon = new((/51,1/),"float",-999.) Nlon = lon ; print(Nlon) ; lon@_FillValue=integertoshort(-999) ;Panel 2 Read in data delim = " " nfields1 = str_fields_count(strs1(0), delim) ; print(nfields1) fieldavg1 = 4 ;avged data R1 = tofloat(str_get_field(strs1,4," ")) R1@_FillValue=integertoshort(-999) print(R) lat1 = tofloat(str_get_field(strs1,2," ")) print(lat) ; Nlat = new((/51,1/),"float",-999.) Nlat1 = lat ; print(Nlat) ; lat1@_FillValue=integertoshort(-999) lon1 = tofloat(str_get_field(strs1,3," ")) print(lon1) ; Nlon = new((/51,1/),"float",-999.) Nlon = lon ; print(Nlon) ; lon@_FillValue=integertoshort(-999) ;Panel 3 Read in data delim = " " nfields2 = str_fields_count(strs2(0), delim) ; print(nfields2) fieldavg2 = 4 ;avged data R2 = tofloat(str_get_field(strs2,4," ")) R2@_FillValue=integertoshort(-999) print(R) lat2 = tofloat(str_get_field(strs,2," ")) print(lat2) ; Nlat = new((/51,1/),"float",-999.) Nlat2 = lat2 ; print(Nlat) ; lat@_FillValue=integertoshort(-999) lon2 = tofloat(str_get_field(strs2,3," ")) print(lon) ; Nlon = new((/51,1/),"float",-999.) Nlon2 = lon2 ; print(Nlon) ; lon2@_FillValue=integertoshort(-999) ;Panel 4 Read in data delim = " " nfields3 = str_fields_count(strs3(0), delim) ; print(nfields) fieldavg3 = 4 ;avged data R3 = tofloat(str_get_field(strs3,4," ")) R3@_FillValue=integertoshort(-999) print(R) lat3 = tofloat(str_get_field(strs3,2," ")) print(lat3) ; Nlat = new((/51,1/),"float",-999.) Nlat3 = lat3 ; print(Nlat) ; lat@_FillValue=integertoshort(-999) lon3 = tofloat(str_get_field(strs3,3," ")) print(lon3) ; Nlon = new((/51,1/),"float",-999.) Nlon3 = lon3 ; print(Nlon) ; lon@_FillValue=integertoshort(-999) ;--------------------------- ; 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/) nlevels = dimsizes(levels) ; colors = read_colormap_file("WhViBlGrYeOrRe") ; setting color markers colors0 = (/"10","30","38","48","56","66","74","94","450","600"/) labels = new(dimsizes(levels)+1,string) ; Labels for legend. ;------------------------------------------------------------- ;Panel 1 info 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 ;------------------------------------------------------------- ;Panel 2 info num_distinct_markers = dimsizes(levels)+1 ; number of distinct markers lat1_new = new((/num_distinct_markers,dimsizes(R1)/),float,-999) lon1_new = new((/num_distinct_markers,dimsizes(R1)/),float,-999) ; print(arr) ; print(lat1_new) ; print(lon1_new) ; ; Group the points according to which range they fall in. At the ; same time, create the label that we will use lat1er in the legend. ; do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(R1.lt.levels(i)) labels(i) = "x < " + levels(i) end if if (i.eq.num_distinct_markers-1) then indexes = ind(R1.ge.max(levels)) labels(i) = "x >= " + max(levels) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(R1.ge.levels(i-1).and.R1.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 lat1/lon1 values and store ; them, so lat1er 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. lat1_new(i,0:npts_range-1) = lat1(indexes) lon1_new(i,0:npts_range-1) = lon1(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do ;------------------------------------------------------------- ;Panel 3 info num_distinct_markers = dimsizes(levels)+1 ; number of distinct markers lat2_new = new((/num_distinct_markers,dimsizes(R2)/),float,-999) lon2_new = new((/num_distinct_markers,dimsizes(R2)/),float,-999) ; print(arr) ; print(lat1_new) ; print(lon1_new) ; ; Group the points according to which range they fall in. At the ; same time, create the label that we will use lat1er in the legend. ; do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(R2.lt.levels(i)) labels(i) = "x < " + levels(i) end if if (i.eq.num_distinct_markers-1) then indexes = ind(R2.ge.max(levels)) labels(i) = "x >= " + max(levels) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(R2.ge.levels(i-1).and.R2.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 lat1/lon1 values and store ; them, so lat1er 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. lat2_new(i,0:npts_range-1) = lat2(indexes) lon2_new(i,0:npts_range-1) = lon2(indexes) end if delete(indexes) ; Necessary b/c "indexes" may be a different ; size next time. end do ;------------------------------------------------------------- ;Panel 4 info num_distinct_markers = dimsizes(levels)+1 ; number of distinct markers lat3_new = new((/num_distinct_markers,dimsizes(R3)/),float,-999) lon3_new = new((/num_distinct_markers,dimsizes(R3)/),float,-999) ; print(arr) ; print(lat1_new) ; print(lon1_new) ; ; Group the points according to which range they fall in. At the ; same time, create the label that we will use lat1er in the legend. ; do i = 0, num_distinct_markers-1 if (i.eq.0) then indexes = ind(R3.lt.levels(i)) labels(i) = "x < " + levels(i) end if if (i.eq.num_distinct_markers-1) then indexes = ind(R3.ge.max(levels)) labels(i) = "x >= " + max(levels) end if if (i.gt.0.and.i.lt.num_distinct_markers-1) then indexes = ind(R3.ge.levels(i-1).and.R3.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 lat1/lon1 values and store ; them, so lat1er 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. lat3_new(i,0:npts_range-1) = lat2(indexes) lon3_new(i,0:npts_range-1) = lon2(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","O3-NYS-DEC-GEOS-Rval4spanel-.png") ; send graphics to PNG file ; gsn_define_colormap(wks,"gui_default") ; define a different color map gsn_define_colormap(wks,"WhViBlGrYeOrRe") ; plot = new(4,graphic) ; Set up some map resources ;---------------Panel 1 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..3" 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 R Winter Values" map0 = gsn_csm_map(wks,mpres) ;---------------Panel 2 mpres1 = True mpres1@gsnDraw = False ; Will draw later mpres1@gsnMaximize = True ; Maximize plot in frame. mpres1@gsnFrame = False ; Don't advance the frame mpres1@pmTickMarkDisplayMode = "Always" ; We want county outlines. mpres1@mpDataBaseVersion = "MediumRes" mpres1@mpDataSetName = "Earth..3" mpres1@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres1@mpOutlineSpecifiers = (/"Land","New York:counties"/) ; Zoom in on NY, U.S. mpres1@mpLimitMode = "LatLon" mpres1@mpMinLatF = 40 mpres1@mpMaxLatF = 46 mpres1@mpMinLonF = -80.0 mpres1@mpMaxLonF = -71.0 ; Turn on tickmarks. mpres1@gsnTickMarksOn = True mpres1@mpFillColors = (/"transparent","transparent", \ "lightgray","transparent"/) ;assign light gray to land masses ; Main title mpres1@tiMainString = "DEC Ozone Site R Spring Values" map1 = gsn_csm_map(wks,mpres1) ;---------------Panel 3 mpres2 = True mpres2@gsnDraw = False ; Will draw later mpres2@gsnMaximize = True ; Maximize plot in frame. mpres2@gsnFrame = False ; Don't advance the frame mpres2@pmTickMarkDisplayMode = "Always" ; We want county outlines. mpres2@mpDataBaseVersion = "MediumRes" mpres2@mpDataSetName = "Earth..3" mpres2@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres2@mpOutlineSpecifiers = (/"Land","New York:counties"/) ; Zoom in on NY, U.S. mpres2@mpLimitMode = "LatLon" mpres2@mpMinLatF = 40 mpres2@mpMaxLatF = 46 mpres2@mpMinLonF = -80.0 mpres2@mpMaxLonF = -71.0 ; Turn on tickmarks. mpres2@gsnTickMarksOn = True mpres2@mpFillColors = (/"transparent","transparent", \ "lightgray","transparent"/) ;assign light gray to land masses ; Main title mpres2@tiMainString = "DEC Ozone Site Summer R Values" map2 = gsn_csm_map(wks,mpres2) ;---------------Panel 4 mpres3 = True mpres3@gsnDraw = False ; Will draw later mpres3@gsnMaximize = True ; Maximize plot in frame. mpres3@gsnFrame = False ; Don't advance the frame mpres3@pmTickMarkDisplayMode = "Always" ; We want county outlines. mpres3@mpDataBaseVersion = "MediumRes" mpres3@mpDataSetName = "Earth..3" mpres3@mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres3@mpOutlineSpecifiers = (/"Land","New York:counties"/) ; Zoom in on NY, U.S. mpres3@mpLimitMode = "LatLon" mpres3@mpMinLatF = 40 mpres3@mpMaxLatF = 46 mpres3@mpMinLonF = -80.0 mpres3@mpMaxLonF = -71.0 ; Turn on tickmarks. mpres3@gsnTickMarksOn = True mpres3@mpFillColors = (/"transparent","transparent", \ "lightgray","transparent"/) ;assign light gray to land masses ; Main title mpres3@tiMainString = "DEC Ozone Site R Fall Values" map3 = gsn_csm_map(wks,mpres3) ; ; 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. ;;-----------Panel 1 mpres = True mpres@gsMarkerIndex = 16 ; Use filled dots for markers. gsres = True gsres@gsMarkerIndex = 16 ; Use filled dots for markers. ;-----------Panel 2 mpres1 = True mpres1@gsMarkerIndex = 16 ; Use filled dots for markers. gsres1 = True gsres1@gsMarkerIndex = 16 ; Use filled dots for markers. ;-----------Panel 3 mpres2 = True mpres2@gsMarkerIndex = 16 ; Use filled dots for markers. gsres2 = True gsres2@gsMarkerIndex = 16 ; Use filled dots for markers. ;-----------Panel 4 mpres3 = True mpres3@gsMarkerIndex = 16 ; Use filled dots for markers. gsres3 = True gsres3@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. ;----------------------------------------------Panel 1 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@gsMarkerThicknessF = 0.7 pmid(i) = gsn_add_polymarker(wks,map0,lon_new(i,:),lat_new(i,:),gsres) end if end do ;----------------------------------------------Panel 2 base_size = 1 pmid1 = new(num_distinct_markers,graphic) do i = 0, num_distinct_markers-1 if (.not.ismissing(lat1_new(i,0))) gsres1@gsMarkerColor = colors(i) gsres1@gsMarkerThicknessF = 0.7*(i+1) gsres1@gsMarkerColor = colors(i) gsres1@gsMarkerThicknessF = 0.7*(i+1) ; gsres1@gsMarkerThicknessF = 0.7 pmid1(i) = gsn_add_polymarker(wks,map1,lon1_new(i,:),lat1_new(i,:),gsres1) end if end do ;----------------------------------------------Panel 3 base_size = 1 pmid2 = new(num_distinct_markers,graphic) do i = 0, num_distinct_markers-1 if (.not.ismissing(lat2_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 pmid2(i) = gsn_add_polymarker(wks,map2,lon2_new(i,:),lat2_new(i,:),gsres) end if end do ;----------------------------------------------Panel 4 base_size = 1 pmid3 = new(num_distinct_markers,graphic) do i = 0, num_distinct_markers-1 if (.not.ismissing(lat3_new(i,0))) gsres@gsMarkerColor = colors3(i) gsres@gsMarkerThicknessF = 0.7*(i+1) gsres@gsMarkerColor = colors3(i) gsres@gsMarkerThicknessF = 0.7*(i+1) ; gsres@gsMarkerThicknessF = 0.7 pmid3(i) = gsn_add_polymarker(wks,map3,lon3_new(i,:),lat3_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,map0,levels,colors0) add_labelbar(wks,map1,levels,colors1) add_labelbar(wks,map2,levels,colors2) add_labelbar(wks,map3,levels,colors3) ; draw(map0) panel_res = True panel_res@gsnPanelRowSpec = True gsn_panel(wks,(/map0,map1,map2,map3/),(/2,2/),panel_res) ; frame(wks) ; Advance the frame. end