Re: drawing data over shapefile still problems

Hi Ugo,

After playing with your script a bit, I'm not really sure you are going to be able to get what you want with the current version of ncl (a future version, yes). That is, if you plot the contour map on top of the shapefile, it will obscure the shapefile. If you plot the shapefile on top of the map, using the solid fill colors in your script, it will obscure the contour plot underneath.

Attached is a plot and modified version of your script. Here, I removed the call to gsn_csm_map(), and moved the shapefile plotting to *after* the gsn_csm_contour_map(). I also made all the solid-fill colors for the shapefile features "transparent". (I also took out the wtaerways and natural shapefiles, as they were taking forever to plot).

Is this something like what you are after? There's no way currently to make colors translucent. Does that help?

Rick

On Apr 29, 2011, at 12:56 PM, ugo merlini wrote:

> Hi,
>
> I'm trying to draw data over shapefile and i have still problems.
>
> here my script
>
> ;----------------------------------------------------------------------
> ; 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
>
> ;--- Open workstation and define colormap.
> wks_type = "png" ; use "newpdf" instead of "pdf" for smaller files
> wks = gsn_open_wks(wks_type,"naturalearth")
>
> ;---Set some resources for the map.
>
> res = True
> res@mpFillOn = False ; No fill will be done.
> res@mpOutlineOn = False ; No outlines will be done.
> res@gsnTickMarksOn = False ; Create contour on the map with lat / lon value
> res@gsnDraw = False ; do not draw the plot
> res@gsnFrame = False ; do not flip the page
> res@gsnMaximize = True ; blow up plot
>
>
> ;---Zoom in on Italy
> res@mpLimitMode = "LatLon"
> res@mpMinLatF = 35.0
> res@mpMaxLatF = 48.0
> res@mpMinLonF = 6.5
> res@mpMaxLonF = 19.0
>
> lnres = True ; resources for polylines
> plres = True ; resources for polygons
>
> map = gsn_csm_map(wks,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)
>
> ;***************************
> ; read in data
> ;***************************
>
> url = "http://nomad1.ncep.noaa.gov:9090/dods/gfs_master/gfs20110427/"
> 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"
>
> res1 = True ; plot mods desired
> res1@cnFillOn = True ; turn on color fill
> res1@cnLinesOn = True ; turn on contour lines
> res1@cnLevelSpacingF = 2 ; contour spacing
> res1@gsnSpreadColors = True ; use full range of color map
> res1@lbLabelStride = 4
> res1@mpOutlineOn = False
> res1@cnSmoothingOn = True
> res1@mpLimitMode="LatLon"
> res1@mpMinLatF = 35.0
> res1@mpMaxLatF = 48.0
> res1@mpMinLonF = 6.5
> res1@mpMaxLonF = 19.0
> res1@gsnDraw = False ; do not draw the plot
> res1@gsnFrame = False ; do not flip the page
> res1@gsnSpreadColors = True ; use full range of color map
> res1@gsnMaximize = True
>
> data_map = gsn_csm_contour_map(wks,TEMP2M,res1)
>
> draw(data_map)
> frame(wks)
> end
>
> any suggestion?
>
> Thanks
>
> Ugo
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Apr 29 16:27:17 2011