;************************************************* ; shapefiles_3.ncl ; ; Concepts illustrated: ; - Reading shapefiles ; - Plotting data from shapefiles ; - Drawing selected data based upon a database query of the shapefile ; - Decreasing the font size of the main title ; - Using data from shapefiles to draw the streams of South America ; - Zooming in on South America on a cylindrical equidistant map ; - Drawing a map using the medium resolution map outlines ; ; Simple example of how to draw selected geometry from a shapefile. ; ; You must download the "HYDRO1k Streams data set for South America ; as a tar file in shapefile format (2.4 MB)" from: ; ; http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/gtopo30/hydro/samerica ; ; Gunzip and untar the file: ; ; gunzip sa_str.tar.gz ; tar -xf sa_str.tar ; ; You then can use NCL V5.1.1 to plot the shapefile data. ; ;************************************************* procedure draw_shapefile_lines(wks,plot) begin f = addfile("JUJUY.shp", "r") ; Open shapefile ; ; Read data off shapefile ; segments = f->segments 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 numFeatures = geomDims(0) ; print(segsDims(0)) ; Number of lines we'll end up drawing. ;************************************************* ; Section to *draw* polylines on existing plot. ;************************************************* plres = True ; resources for polylines plres@gsLineColor = "blue" plres@gsLineThicknessF = 2.0 lon = f->x lat = f->y segNum = 0 ; Counter for adding polylines ;---Loop through each segment and attach to plot 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 gsn_polyline(wks, plot, lon(startPT:endPT), \ lat(startPT:endPT), plres) segNum = segNum + 1 end do end do end