load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" ;---------------------------------------------------------------------- ; This procedure attaches a box to a map. ; ; It gets the edges of the box using gc_latlon to calculate the ; great circle between two lat/lon points. ; ; This ensures that your box edges will be straight when you attach ; them to the map. ;---------------------------------------------------------------------- function add_box(wks,map,lftlat,lftlon,rgtlat,rgtlon) local lnres, npts, i, lat_begend, lon_begend, dist begin lnres = True lnres@gsLineColor = "red" lnres@gsLineThicknessF = 2.0 npts = 20 ; Number of points along each box edge. ; You could make this different for each ; edge, if you want. ; ; Define the coordinates for the start, end of the four sides ; of each box. Put them in a big array so it's easier to ; loop across the points later. ; ; bottom right ; top left lat_begend = (/ (/lftlat,lftlat/), (/lftlat,rgtlat/), \ (/rgtlat,rgtlat/), (/rgtlat,lftlat/)/) lon_begend = (/ (/lftlon,rgtlon/), (/rgtlon,rgtlon/), \ (/rgtlon,lftlon/), (/lftlon,lftlon/)/) ;---Define array to hold box. latbox = new(4*npts,float) lonbox = new(4*npts,float) ;---Loop across the four edges and calculate the points along each edge do i=0,3 ibeg = i*npts iend = ibeg+npts-1 dist = gc_latlon(lat_begend(i,0),lon_begend(i,0), \ lat_begend(i,1),lon_begend(i,1),npts,2) latbox(ibeg:iend) = dist@gclat lonbox(ibeg:iend) = dist@gclon end do ;---This is so we can cross the 0 longitude correctly. lonbox = where(lonbox.gt.180,lonbox-360,lonbox) ;---Attach box to map and return dum = gsn_add_polyline(wks, map, lonbox, latbox, lnres) return(dum) end ;---------------------------------------------------------------------- ; Main code ;---------------------------------------------------------------------- begin wks = gsn_open_wks("ps", "boxes") ; Open PostScript file. res = True res@gsnDraw = False res@gsnFrame = False res@gsnMaximize = True ; Add map resources res@mpOutlineDrawOrder = "PostDraw" ; Draw map outlines last res@mpGridAndLimbOn = False ; Turn off lat/lon lines res@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks res@mpOutlineDrawOrder = "PostDraw" ; Draw map outlines last res@mpGridAndLimbOn = False ; Turn off lat/lon lines res@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks res@mpLimitMode = "Corners" ; Portion of map to zoom res@mpLimitMode = "Corners" ; Portion of map to zoom res@mpLeftCornerLatF = 25. res@mpLeftCornerLonF = -14. res@mpRightCornerLatF = 55. res@mpRightCornerLonF = 48. res@tiMainString = "Box added using gsn_add_polyline" map1 = gsn_csm_map(wks, res) ;---Draw a box using gsn_add_polyline. See how the edges are curved. lnres = True lnres@gsLineColor = "Blue" lnres@gsLineThicknessF = 2.0 lftlat = 27.6816 rgtlat = 46.6198 lftlon = -12.9163 rgtlon = 44.2297 lonbox = (/ lftlon, rgtlon, rgtlon, lftlon, lftlon/) latbox = (/ lftlat, lftlat, rgtlat, rgtlat, lftlat/) dum1 = gsn_add_polyline(wks, map1, lonbox, latbox, lnres) draw(map1) frame(wks) ;---Create new map plot with a new title. res@tiMainString = "Box added using gc_latlon" map2 = gsn_csm_map(wks, res) ;---Calculate box edges using gc_latlon. dum2 = add_box(wks,map2,lftlat,lftlon,rgtlat,rgtlon) draw(map2) frame(wks) ;---Create new map plot with a new title. res@tiMainString = "Box added using gsn_blank_plot" map3 = gsn_csm_map(wks, res) ;---Get viewport coordiates. getvalues map3 "vpXF" : vpx "vpYF" : vpy "vpWidthF" : vpw "vpHeightF" : vph end getvalues xlndc = new(1,float) ylndc = new(1,float) xrndc = new(1,float) yrndc = new(1,float) datatondc(map3, lftlon, lftlat, xlndc, ylndc) datatondc(map3, rgtlon, rgtlat, xrndc, yrndc) ; ; Create a blank plot that's in same physical location as map3. ; Make sure the axis limits are in the same NDC coordinates as the ; location of the map. ; bres = True bres@trXMinF = vpx bres@trYMinF = vpy-vph bres@trXMaxF = vpx+vpw bres@trYMaxF = vpy bres@pmTickMarkDisplayMode = "Never" ;---This will put blank plot and map in same position. bres@tfDoNDCOverlay = True ;---Create the blank plot with no tickmarks. blank = gsn_blank_plot(wks,bres) ;---Add line to the blank plot. xbox = (/ xlndc, xrndc, xrndc, xlndc, xlndc/) ybox = (/ ylndc, ylndc, yrndc, yrndc, ylndc/) lnres@gsLineColor = "Purple" dum2 = gsn_add_polyline(wks, blank, xbox, ybox, lnres) ;--Overlay blank plot on map, and you will see the line. overlay(map3,blank) draw(map3) frame(wks) end