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" undef("create_map") function create_map(wks,lat,lon) local gres, lnres, f, map, segments, geometry, segsDims, geomDims, \ geom_segIndex, geom_numSegs, segs_xyzIndex, segs_numPnts, numFeatures begin ; Adjust limits to plot minlat = min(lat)-2 maxlat = max(lat)+4 minlon = min(lon)-3 maxlon = max(lon)+1 ;---Open shapefile and read data f = addfile("COL_adm1.shp","r") segments = f->segments geometry = f->geometry segsDims = dimsizes(segments) geomDims = dimsizes(geometry) geom_segIndex = f@geom_segIndex geom_numSegs = f@geom_numSegs segs_xyzIndex = f@segs_xyzIndex segs_numPnts = f@segs_numPnts numFeatures = geomDims(0) lon1 = f->x lat1 = f->y ;---Set some map resources gres = True gres@gsnDraw = False gres@gsnFrame = False gres@mpOutlineOn = False gres@mpFillOn = False gres@mpPerimOn = True gres@pmTickMarkDisplayMode = "Always" gres@mpProjection = "Mercator" gres@mpLimitMode = "LatLon" gres@mpMinLatF = min((/min(lat1),minlat/)) gres@mpMaxLatF = max((/max(lat1),maxlat/)) gres@mpMinLonF = min((/min(lon1),minlon/)) gres@mpMaxLonF = max((/max(lon1),maxlon/)) ; print("min/max lat = " + minlat + "/" + maxlat) ; print("min/max lon = " + minlon + "/" + maxlon) ; print("min/max lat1 = " + min(lat1) + "/" + max(lat1)) ; print("min/max lon1 = " + min(lon1) + "/" + max(lon1)) ;---Create the map. map = gsn_csm_map(wks,gres) ;---Set some line resources. lnres = True lnres@gsLineThicknessF = 1.0 lnres@gsLineColor = "black" ;---Attach the shapefile polylines. do i=0, numFeatures-1 startSegment = geometry(i, geom_segIndex) numSegments = geometry(i, geom_numSegs) do seg=startSegment, startSegment+numSegments-1 startPT = segments(seg, segs_xyzIndex) endPT = startPT + segments(seg, segs_numPnts) - 1 dumstr = unique_string("primitive") map@$dumstr$ = gsn_add_polyline(wks, map, lon1(startPT:endPT), \ lat1(startPT:endPT), lnres) end do end do ;---Return map to main code. return(map) end ;---------------------------------------------------------------------- ; Main code ;---------------------------------------------------------------------- begin wks = gsn_open_wks("ps","com_rf") gsn_define_colormap(wks,"BlWhRe") gsn_reverse_colormap(wks) fname = "pr-ene.txt" lines = asciiread(fname,-1,"string") delim = str_get_tab() lat = stringtofloat(str_get_field(lines(1:),2,delim)) lon = stringtofloat(str_get_field(lines(1:),3,delim)) idx = tointeger(str_get_field(lines(1:),7,delim)) map = create_map(wks,lat,lon) res = True res@gsnDraw = False res@gsnFrame = False res@gsnSpreadColors = True res@gsnSpreadColorStart = 24 res@gsnSpreadColorEnd = -26 res@cnLineLabelPlacementMode = "Constant" res@cnLineDashSegLenF = 0.00003 res@cnLevelSelectionMode = "AutomaticLevels" res@cnFillOn = True res@cnLinesOn = False res@cnLineLabelsOn = False res@cnLevelFlags = new(139,"string") res@cnLevelFlags(:) = "NoLine" res@cnLevelFlags(0::20) = "LineAndLabel" res@cnSmoothingOn = True res@lbBoxLinesOn = False res@tiMainString = "RAINFALL - JAN 2011" res@tiMainFontHeightF = 0.02 res@lbLabelFontHeightF = 0.01 res@lbTitleOn = True res@lbTitleString = "Rainfall Probability" res@lbTitleFontHeightF = 0.01 res@sfXArray = lon res@sfYArray = lat ;---Create contours and overlay them on map contour = gsn_csm_contour(wks,idx,res) overlay(map,contour) ;---Maximize new output in PostScript frame. pres = True maximize_output(wks,pres) end