;************************************************* ; shapefiles_1.ncl ; ; Concepts illustrated: ; - Reading shapefiles ; - Plotting data from shapefiles ; - Using shapefile data to plot unemployment percentages in the U.S. ; - Drawing a custom labelbar on a map ; - Drawing filled polygons over a Lambert Conformal plot ; - Drawing the US with a Lambert Conformal projection ; - Zooming in on a particular area on a Lambert Conformal map ; - Centering the labels under the labelbar boxes ; ;************************************************* ; ; Simple example of how to draw selected geometry from a shapefile, ; based upon properties of an associated non-spatial variable. ; ; This example color-fills the states based upon "percent unemployment", ; which is calculated from several of the non-spatial variables in the ; file. ; ; "states.shp" is from the National Atlas (http://www.nationalatlas.gov/) ; ; You must also have the files "states.dbf" and "states.shx" for this ; example to run. ;************************************************* ; These files are loaded by default in NCL V6.2.0 and newer ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" begin wks = gsn_open_wks("eps","shapefiles") ; send graphics to PNG file res = True res@gsnDraw = False ; don't draw yet res@gsnFrame = False ; don't advance frame yet res@gsnMaximize = True ; maximize plot in frame ;res@mpProjection = "LambertConformal" ; choose projection res@mpLambertParallel1F = 33 ; first parallel res@mpLambertParallel2F = 45 ; second parallel res@mpLambertMeridianF = -98 ; meridian res@mpLimitMode = "Corners" ; corner method of zoom res@mpLeftCornerLatF = 29 ; left corner res@mpLeftCornerLonF = -8 ; left corner res@mpRightCornerLatF = 48 ; right corner res@mpRightCornerLonF = 45 ; right corner res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks res@tiMainString = "GSAs" plot = gsn_csm_map(wks,res) ; Create map, but don't draw it yet. ;************************************************* ; Section to add polygons to map. ;************************************************* dir = ncargpath("data") + "/shp/" ; dir = ncargpath("/home/netapp-clima/scratch/asalhi/mask/") + "/shp/" ;f = addfile(dir + "states.shp", "r") ; Open shapefile f = addfile("/home/netapp-clima/scratch/asalhi/mask/GSAs_simplified.shp", "r") ; Open shapefile ; ; Read data off shapefile ; segments = f->segments ;segments = f->CENTER_X,N,24,15 ;geometry = f->CENTER_yN2415 geometry = f->geometry segsDims = dimsizes(segments) geomDims = dimsizes(geometry) ; ; Read global attributes ; geom_segIndex = f@geom_segIndex geom_numSegs = f@geom_numSegs segs_xyzIndex = f@segs_xyzIndex segs_numPnts = f@segs_numPnts lines = new(segsDims(0),graphic) ; Array to hold polygons numFeatures = geomDims(0) ;unemp = f->UNEMPLOY / f->PERSONS ; GSA = f->CENTER_X /f->CENTER_Y GSA = f->SMU_CODE ; GSA = f->SECT_CODC80 lon = f->x lat = f->y ;lon = f->CENTER_X ; lat = f->CENTER_Y plres = True ; resources for polylines plres@gsEdgesOn = True ; draw border around polygons plres@gsEdgeColor = "black" colors = (/"blue","green","yellow","red","MediumTurquoise","MediumSpringGreen","HotPink","violet","magenta","honeydew4","LightCyan1","khaki3","salmon3","orange3","maroon3","VioletRed4","gray7","papaya whip","lemon chiffon","LavenderBlush","MistyRose","MidnightBlue","aquamarine","SandyBrown","coral","LightPink","ivory2","chartreuse2","goldenrod3","chocolate2","firebrick3","linen"/) segNum = 0 do i=0, numFeatures-1 ; color assignment (probably a better way to do this?) if (GSA(i).ge.0.01 .and. GSA(i).lt.0.02) then plres@gsFillColor = colors(0) end if if (GSA(i).ge.0.02 .and. GSA(i).lt.0.03) then plres@gsFillColor = colors(1) end if if (GSA(i).ge.0.03 .and. GSA(i).lt.0.04) then plres@gsFillColor = colors(2) end if if (GSA(i).ge.0.04 .and. GSA(i).lt.0.05) then plres@gsFillColor = colors(4) end if if (GSA(i).ge.0.05 .and. GSA(i).lt.0.06) then plres@gsFillColor = colors(5) end if if (GSA(i).ge.0.06 .and. GSA(i).lt.0.07) then plres@gsFillColor = colors(6) end if if (GSA(i).ge.0.07 .and. GSA(i).lt.0.08) then plres@gsFillColor = colors(7) end if if (GSA(i).ge.0.08 .and. GSA(i).lt.0.09) then plres@gsFillColor = colors(8) end if if (GSA(i).ge.0.09 .and. GSA(i).lt.0.010) then plres@gsFillColor = colors(9) end if if (GSA(i).ge.0.010 .and. GSA(i).lt.0.11) then plres@gsFillColor = colors(10) end if if (GSA(i).ge.0.011 .and. GSA(i).lt.0.0111) then plres@gsFillColor = colors(11) end if if (GSA(i).ge.0.0111 .and. GSA(i).lt.0.0112) then plres@gsFillColor = colors(111) end if if (GSA(i).ge.0.0112 .and. GSA(i).lt.0.012) then plres@gsFillColor = colors(112) end if if (GSA(i).ge.0.012 .and. GSA(i).lt.0.013) then plres@gsFillColor = colors(12) end if if (GSA(i).ge.0.013 .and. GSA(i).lt.0.014) then plres@gsFillColor = colors(13) end if if (GSA(i).ge.0.014 .and. GSA(i).lt.0.015) then plres@gsFillColor = colors(14) end if if (GSA(i).ge.0.015 .and. GSA(i).lt.0.016) then plres@gsFillColor = colors(15) end if if (GSA(i).ge.0.016 .and. GSA(i).lt.0.017) then plres@gsFillColor = colors(16) end if if (GSA(i).ge.0.017 .and. GSA(i).lt.0.018) then plres@gsFillColor = colors(17) end if if (GSA(i).ge.0.018 .and. GSA(i).lt.0.019) then plres@gsFillColor = colors(18) end if if (GSA(i).ge.0.019 .and. GSA(i).lt.0.020) then plres@gsFillColor = colors(19) end if if (GSA(i).ge.0.020 .and. GSA(i).lt.0.021) then plres@gsFillColor = colors(20) end if if (GSA(i).ge.0.021 .and. GSA(i).lt.0.022) then plres@gsFillColor = colors(21) end if if (GSA(i).ge.0.022 .and. GSA(i).lt.0.023) then plres@gsFillColor = colors(22) end if if (GSA(i).ge.0.023 .and. GSA(i).lt.0.024) then plres@gsFillColor = colors(23) end if if (GSA(i).ge.0.024 .and. GSA(i).lt.0.025) then plres@gsFillColor = colors(24) end if if (GSA(i).ge.0.025 .and. GSA(i).lt.0.026) then plres@gsFillColor = colors(25) end if if (GSA(i).ge.0.026 .and. GSA(i).lt.0.027) then plres@gsFillColor = colors(26) end if if (GSA(i).ge.0.027 .and. GSA(i).lt.0.028) then plres@gsFillColor = colors(27) end if if (GSA(i).ge.0.028 .and. GSA(i).lt.0.029) then plres@gsFillColor = colors(28) end if if (GSA(i).ge.0.029 .and. GSA(i).lt.0.030) then plres@gsFillColor = colors(29) end if if (GSA(i).ge.0.030) then plres@gsFillColor = colors(30) end if 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 lines(segNum) = gsn_add_polygon(wks, plot, lon(startPT:endPT), \ lat(startPT:endPT), plres) segNum = segNum + 1 end do end do ; Make a legend... labels = (/ "1", "2", "3", "4","","5","6","7","8","9","10","11","11.1","11.2","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"/) lres = True lres@vpWidthF = 0.50 ; location lres@vpHeightF = 0.05 ; " " lres@lbPerimOn = False ; Turn off perimeter. lres@lbOrientation = "Horizontal" ; Default is vertical. lres@lbLabelAlignment = "BoxCenters" ; Default is "BoxCenters". lres@lbFillColors = colors lres@lbMonoFillPattern = True ; Fill them all solid. lres@lbLabelFontHeightF = 0.012 ; label font height ;lres@lbTitleString = "percent" ; title lres@lbTitlePosition = "Bottom" ; location of title lres@lbTitleFontHeightF = 0.01 ; title font height gsn_labelbar_ndc (wks,32,labels,0.23,0.15,lres) ; ; Maximize output in frame. This will draw everything and advance ; the frame. ; maximize_output(wks,False) end