;---------------------------------------------------------------------- ; shapefiles_6.ncl ; ; Concepts illustrated: ; - Reading shapefiles ; - Plotting data from shapefiles ; - Using data from shapefiles to draw the cantons of Switzerland ; - Zooming in on Switzerland on a cylindrical equidistant map ; - Creating a color map using named colors ; - Attaching lots of text strings to a map ;---------------------------------------------------------------------- ; This example shows how to read geographic data from a shapefile ; and plot it on a map created by NCL. ;---------------------------------------------------------------------- ; This particular example plots data for Switzerland. ;---------------------------------------------------------------------- ; Got the "switzerland" shapefiles directory from: ; ; http://download.geofabrik.de/osm/europe/ ; ; Got the "CHE_adm" shapefiles directory from: ; ; http://www.gadm.org/country ;---------------------------------------------------------------------- load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" begin ;--- Open workstation and define colormap. wks = gsn_open_wks("pdf","JUJUYshp") gsn_define_colormap(wks,(/"white","black","tan","LightBlue","brown","yellow"/)) ;---Set some resources for the map. res = True res@gsnMaximize = True res@mpFillOn = False res@mpOutlineBoundarySets = "AllBoundaries" res@gsnTickMarksOn = False ;---Zoom in on Swizterland res@mpLimitMode = "LatLon" res@mpMinLatF = -25 res@mpMaxLatF = -23 res@mpMinLonF = -66 res@mpMaxLonF = -63 ;---Draw the crude map of Switzerland. map = gsn_csm_map(wks,res) ;---Read shapefile with "places" information. f = addfile("JUJUY-Points.shp","r") print("Reading 'JUJUY-Points.shp'") prims = True names = f->NAME lon = f->x lat = f->y num_places = dimsizes(names) ;---Set up text resources to label random areas of interest. txres = True txres@txFontHeightF = 0.008 txres@txJust = "CenterLeft" ;---Set up marker resources to mark random areas of interest. mkres = True mkres@gsMarkerIndex = 17 mkres@gsMarkerColor = "green" mkres@gsMarkerSizeF = 30. ;---Areas we want to put a label and a marker. names_of_interest = (/"Aeroclub","AEROP Salta","Gral.Güemes","Perico","Lib.Gral.San Martin",\ "Campo Quijano","Rosario de Lerma","El Carril","San Pedro","Cerrillos",\ "Chicoana","Cnel.Moldes","Salta","La Viña","S.S.Jujuy",\ "El Carmen","Monterrico"/) ; ; Loop through each of the "name" places in the shapefile, and see if ; it's on our list of ones we want to put a label and marker for. ; ; Also, make sure we haven't already encountered this place. The ; shapefile seems to contain duplicate entries (Zurich has two entries, ; for example). ; ; Just like with gsn_add_polyline above, if you find this looping ; code to be too slow, you can use gsn_polymarker and gsn_text instead. ; names_of_interest@_FillValue = "missing" do i=0,num_places-1 if(any(names(i).eq.names_of_interest)) then ii = ind(names_of_interest.eq.names(i)) names_of_interest(ii) = "missing" dumstr = unique_string("primitive") ;---Attach the marker to the map. print (lon(i)) print (lat(i)) prims@$dumstr$ = gsn_add_polymarker(wks, map, lon(i), lat(i), mkres) dumstr = unique_string("primitive") ; ; NCL doesn't deal with "ü" unfortunately! This code should be done ; outside this loop. ; if(names(i).eq."Gral.Güemes") then names(i) = "Gral.Guemes" end if ;---Attach the text string to the map. prims@$dumstr$ = gsn_add_text(wks, map, " " + names(i), \ lon(i), lat(i), txres) delete(ii) end if end do ;---We're done adding primitives, draw everything and advance frame. draw(map) frame(wks) end