Re: ?m able to get the data but merge in a shapefile map

Hi Ugo,
If I understand things correctly, you have two problems:
You are getting two pages of output.

To address this problem, you need to tell NCL not to draw the plot or
flip the page when you create the contour plot. To accomplish this, add
the following to your res resource list:

res@gsnDraw = False ; do not draw the plot
res@gsnFrame = False ; do not flip the page

The second issue is as follows: On the first page is the contour plot,
and on the second page is the shapefile polygons overlaid on the contour
plot. The polygons are drawn over the contour plot. You would like the
polygons to be drawn before the contours are drawn.

This is a draw order issue. You need to tell NCL when to draw the
contours, and when to draw the polygons. See this page here:
http://www.ncl.ucar.edu/Applications/draworder.shtml

Specifically, you will need to set tfPolyDrawOrder and you may need to
set cnLineDrawOrder.

To tell NCL to draw the polygons last, add this line to your polygon
resource list:
plres@tfPolyDrawOrder = "PostDraw"

You may have to tell NCL when to draw the contours. To do that set this
in your contour resource list:
res@cnLineDrawOrder = "Draw"

Hope that helps. If it doesn't, let the group know..
Adam

On 04/26/2011 04:26 PM, ugo merlini wrote:
> Hi,
>
> I'm able to grab the data from a file and built a contour map, I'm able
> to create a shapefile map but not merge them....
>
> the script included generate two files one with only the data no map and
> a second one like the first but with the map above (where there is the
> map the data are not displayed)
>
> Where I'm wrong? Suggested reading?
>
> Regards
>
> Ugo
>
>
> ;----------------------------------------------------------------------
> ; This example shows how to read geographic data
> ; from Natural Earth shapefiles
> ; and plot them as polylines and polygons.
> ;----------------------------------------------------------------------
> ; This particular example plots data for Switzerland.
> ;----------------------------------------------------------------------
> ; Download the shapefiles from http://www.naturalearthdata.com/
> ; Unzip to a directory
> ;----------------------------------------------------------------------
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>
> begin
> ;***************************
> ; read in data
> ;***************************
>
> url = "http://nomad1.ncep.noaa.gov:9090/dods/gfs_master/gfs20110420/"
> filename = url + "gfs_master_00z"
>
> exists = isfilepresent(filename)
> if(.not.exists) then
> print("OPeNDAP isfilepresent test unsuccessful.")
> print("Either file doesn't exist, or NCL does not have OPeNDAP
> capabilities on this system")
> else
> print("OPeNDAP isfilepresent test successful.")
> a = addfile(filename,"r")
> vars = getfilevarnames(a)
> end if
>
> TEMP2M = a->tmp2m(1,:,:) ; temperature
> TEMP2M = TEMP2M - 273.15
> TEMP2M@units = "(C)"
> TEMP2M@long_name = "Temperatura a 2 metri dal suolo"
>
> res = True ; plot mods desired
> res@cnFillOn = True ; turn on color fill
> res@cnLinesOn = False ; turn of contour lines
> res@cnLevelSpacingF = 0.5 ; contour spacing
> res@gsnSpreadColors = True ; use full range of color map
> res@lbLabelStride = 4
> res@mpOutlineOn = False
>
> res@mpLimitMode="LatLon"
> res@mpMinLatF = 35.0
> res@mpMaxLatF = 48.0
> res@mpMinLonF = 6.5
> res@mpMaxLonF = 19.0
> res@gsnSpreadColors = True ; use full range of color map
> res@gsnMaximize = True
>
>
> ;--- Open workstation and define colormap.
> wks_type = "png" ; use "newpdf" instead of "pdf" for smaller files
> wks = gsn_open_wks(wks_type,"naturalearth_con")
> gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap
> lnres = True ; resources for polylines
> plres = True ; resources for polygons
>
> map = gsn_csm_contour_map_ce(wks,TEMP2M,res)
>
> ;---Shapefiles
> fnames = "shapefile/italy/" +
> (/"ITA_adm0","ITA_adm1","ITA_adm2","waterways","natural"/) + ".shp"
> linecolors = (/"Black","Black","Black","Blue","Blue"/)
> fillcolors = (/"White","White","White","Blue","Blue"/)
> thick = (/0,0,0,0,0/)
>
> ;---Loop through files that we want to read geographic information from.
> prims = True
> lines = True
> do n=0,dimsizes(fnames)-1
> ;---Open the shapefile.
> f = addfile(fnames(n),"r")
>
> ;---Read data off the 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
> geometry_type = f@geometry_type
>
> numFeatures = geomDims(0)
>
> lon = f->x
> lat = f->y
>
> ; put if statement outside the loop
> if (geometry_type.eq."polygon") then
> plres@gsFillColor = fillcolors(n)
> plres@gsEdgesOn = True ; draw border around polygons
> plres@gsEdgeColor = linecolors(n)
> plres@gsEdgeThicknessF = thick(n)
>
> ;---Section to draw polygons on map.
> 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
> ;---This call adds the polygon.
> dumstr = unique_string("lines")
> lines@$dumstr$ = gsn_add_polygon(wks, map , lon(startPT:endPT),
> lat(startPT:endPT), plres)
> end do
> end do
> else
> lnres@gsLineThicknessF = thick(n)
> lnres@gsLineColor = linecolors(n)
> ;---Section to draw polylines on map.
> 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
> ;---This call adds the line segment.
> dumstr = unique_string("primitive")
> prims@$dumstr$ = gsn_add_polyline(wks, map, lon(startPT:endPT),
> lat(startPT:endPT), lnres)
> end do
> end do
> end if
>
> ;---Clean up before we read in same variables again.
> delete(lat)
> delete(lon)
> delete(segments)
> delete(geometry)
> delete(segsDims)
> delete(geomDims)
> end do
>
> ;---We're done adding primitives, draw everything and advance frame.
>
> draw(map)
> frame(wks)
> end
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
__________________________________________________
Adam Phillips 
asphilli@ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
Climate and Global Dynamics Division         fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000    http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk