drawing data over shapefile still problems

From: ugo merlini <ugomerlini_at_nyahnyahspammersnyahnyah>
Date: Fri Apr 29 2011 - 12:56:42 MDT

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
Received on Fri Apr 29 12:57:54 2011

This archive was generated by hypermail 2.1.8 : Tue May 03 2011 - 14:47:34 MDT