;---------------------------------------------------------------------- ; This series of examples shows how to add shapefile outlines to ; an NCL/contour map. ; ; The 1' topo grid can be downloaded from: ; ; You can download the Poland (and most other countries) ; shapefile data from: ; ; http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/cell_registered/netcdf/ ; ; ; http://gadm.org/country/ ; ;---------------------------------------------------------------------- ;---------------------------------------------------------------------- ; This function reads in a color table and levels provided by: ; http://netgis.geo.uw.edu.pl/srtm/wizualizacja.shtml ; ; Since both levels and a color map are returned by this function, ; we use the special List type to return both variables. ;---------------------------------------------------------------------- undef("earth_model_color_map") function earth_model_color_map() local lines, nlines begin lines = asciiread("Poland_map_and_levels.txt",-1,"string") nlines = dimsizes(lines) levels = tofloat(str_get_field(lines(1:),1,":")) cmap = new((/nlines-1,3/),float) ;---Skip first level and color, because the color is the same cmap(::,0) = tofloat(str_get_field(lines(1:),2,":"))/255. cmap(::,1) = tofloat(str_get_field(lines(1:),3,":"))/255. cmap(::,2) = tofloat(str_get_field(lines(1:),4,":"))/255. return([/cmap,levels/]) end ;---------------------------------------------------------------------- ; Main driver code ;---------------------------------------------------------------------- begin shp_dir = "POL_adm/" shp_filename = "POL_adm1.shp" ; use POL_adm2.shp for more outlines ; ; Get colormap and levels to use. This function returns two variables ; as an NCL "List" object. ; ret = earth_model_color_map() cmap = ret[0] levels = ret[1] ;---Map area to zoom in on minlat = 48.5 maxlat = 55 minlon = 14 maxlon = 24.5 ;---Read subset of data from a 1' file. topo_file = "ETOPO1_Bed_c_gmt4.grd.nc" a = addfile(topo_file,"r") elev = a->z({minlat:maxlat},{minlon:maxlon}) ;---Start the graphics wks = gsn_open_wks("png","shapefile1d") res = True res@gsnDraw = False ; Turn off draw res@gsnFrame = False ; Turn off frame res@gsnMaximize = True ; maximize plot in frame res@cnFillOn = True ; turn on contour fill res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off line labels res@cnInfoLabelOn = False ; turn off info label res@cnFillPalette = cmap ; change color map res@cnLevelSelectionMode = "ExplicitLevels" res@cnLevels = levels res@gsnAddCyclic = False ;---Zoom in on map res@mpMinLatF = minlat res@mpMaxLatF = maxlat res@mpMinLonF = minlon res@mpMaxLonF = maxlon res@mpCenterLonF = (res@mpMinLonF + res@mpMaxLonF) / 2. res@mpOutlineOn = False ; turn off NCL's map outlines res@mpFillOn = False ; turn off map fill res@tiMainString = topo_file res@pmTitleZone = 4 ; Move main title down a little res@gsnLeftString = "elevation" res@gsnRightString = "m" res@pmTickMarkDisplayMode = "Always" ; turn on "nice" tickmark labels plot = gsn_csm_contour_map(wks,elev,res) ;---Attach outlines from shapefile. lnres = True lnres@gsLineThicknessF = 2.0 ; 1.0 is the default dum = gsn_add_shapefile_polylines(wks,plot,shp_dir+shp_filename,lnres) ;---Draw and advance frame draw(plot) frame(wks) end