blank contour over grenwich, why? and questions about limit opendap latitude and longitute search and projection plot

hi

I create this 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,"JAL")
  gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap

  res = True ; plot mods desired
  res@gsnMaximize = True
  ;res@cnFillOn = False ; turn on color fill
  res@gsnTickMarksOn = False
  res@mpLimitMode="LatLon"
  res@mpMinLatF = 27.0
  res@mpMaxLatF = 90.0
  res@mpMinLonF = -74
  res@mpMaxLonF = 52
  res@gsnDraw = False ; do not draw the plot
  res@gsnFrame = False ; do not flip the page
  
   map = gsn_csm_map(wks,res)

;---Shapefiles
  fnames = "/mnt/internetserver/map/shapefile/europa/aaa_full/" + (/"europa"/) + ".shp" ; Files to be open
  linecolors = (/"Black"/) ; color definition for each shapefile
  fillcolors = (/"Transparent"/) ; Fill color definition for each shapefile
  thicks = (/2/) ; Thickness for each shapefile
  lnres = True ; resources for polylines
  plres = True

;---Loop through files that we want to read geographic information from.
  prims = True
  lines = True
  do n=0,dimsizes(fnames)-1 ; Loop through files that we want to read geographic information from.
   f = addfile(fnames(n),"r") ; Open the shapefile number n.
   segments = f->segments ; Read data off the shapefile
   geometry = f->geometry
   segsDims = dimsizes(segments)
   geomDims = dimsizes(geometry)
   geom_segIndex = f@geom_segIndex ; Read global attributes
   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
   if (geometry_type.eq."polygon") then ; Put if statement outside the loop
      plres@gsFillColor = fillcolors(n)
      plres@gsEdgesOn = True ; Draw border around polygons
      plres@gsEdgeColor = linecolors(n)
      plres@gsEdgeThicknessF = thicks(n)
      do i=0, numFeatures-1 ; Section to draw polygons on map.
        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
          dumstr = unique_string("lines") ; This call adds the polygon.
          map@$dumstr$ = gsn_add_polygon(wks, map , lon(startPT:endPT), lat(startPT:endPT), plres)
        end do
      end do
   else
      lnres@gsLineThicknessF = thicks(n)
      lnres@gsLineColor = linecolors(n)
      do i=0, numFeatures-1 ; Section to draw polylines on map.
        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
          dumstr = unique_string("primitive") ; This call adds the line segment.
          map@$dumstr$ = gsn_add_polyline(wks, map, lon(startPT:endPT), lat(startPT:endPT), lnres)
        end do
      end do
   end if
   delete(lat) ; Clean up before we read in same variables again.
   delete(lon)
   delete(segments)
   delete(geometry)
   delete(segsDims)
   delete(geomDims)
  end do

;***************************
; read in data
;***************************

  url = "http://nomad1.ncep.noaa.gov:9090/dods/gfs_master/gfs20110815/"
  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(:,:,:) ; 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
  res1@cnLinesOn = True ; turn contour lines
  res1@gsnSpreadColors = True ; use entire color map
  res1@lbLabelAutoStride = True ; let NCL determine label spacing

  res1@gsnDraw = False ; do not draw the plot
  res1@gsnFrame = False ; do not flip the page
  res1@cnLineDrawOrder = "PostDraw"
  res1@gsnSpreadColors = True ; use full range of color map
  res1@gsnMaximize = True
  i = 0
   do while(i.le.1)

  plot = gsn_csm_contour(wks,TEMP2M(i,:,:),res1)
  overlay(map,plot)
  
;---We're done adding primitives, draw everything and advance frame.
  
  draw(map)
  frame(wks)
  i=i+1
 end do

end

I got the picture like the linked where no contour in draw over greenwich.

http://i103.photobucket.com/albums/m136/ugo73/JAL000001.png

now i took 35 second to parse the shape file but the reading the opendap file is still slow is there a way to limit the longitute latitude reading in it?

if I change the projetion in res eresource it will change automatically in res1?

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Aug 16 07:45:26 2011