(no subject)

i made my first script which use a shapefile as source but it is very slow to
produce the output

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

;*******************

;* Build shapefile *

;*******************

;* Procedure definition *

undef("attach_shapefile_outlines")

procedure attach_shapefile_outlines(wks,map)

local lnres, plres, f, segments, geometry, segsDims, geomDims, fnames,
fillcolors, thicks

;* Procedure code *

begin

;* Shapefile parameters definition

  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
                       ; resources for polygons

;* Shapefiles reading

  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

end

begin

;*************************************

; Open workstation and define colormap

;*************************************

  wks_type = "png"

  wks = gsn_open_wks(wks_type,"pro")

  gsn_define_colormap(wks,"wh-bl-gr-ye-re") ; choose colormap

;******************

;* Read data file *

;******************

  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"

  res = True ; plot mods desired

  res@cnFillOn = True ; turn on/off color between
contour levels

  res@cnLevelSpacingF = 1.0 ; contour spacing

  res@cnLinesOn = True ; turn on/off contour lines
between contour levels

  res@cnLineColor = "White" ; Set color of the line

  res@cnLineThicknessF = 2.0 ; Set the thickness of the line

  res@mpShapeMode = "FreeAspect"

  res@vpWidthF = 0.8

  res@vpHeightF = 0.8

  res@mpLimitMode="Latlon"

  res@mpMinLatF = 27.0

  res@mpMaxLatF = 90.0

  res@mpMinLonF = -74.0

  res@mpMaxLonF = 52.0

  res@mpOutlineOn = False ;turn On/Off standard NCL map

  ;res@gsnMaximize = True

  res@gsnDraw = False ; do not draw the plot

  res@gsnFrame = False ; do not flip the page

  res@lbLabelStride = 5 ; Label spacing

  ;res@gsnPaperOrientation = "auto"

  res@pmTickMarkDisplayMode = "NoCreate"

 i = 0

  do while(i.le.50)

  wks = gsn_open_wks(wks_type,"pro"+i)

  gsn_define_colormap(wks,"wh-bl-gr-ye-re") ; choose colormap

  plot = gsn_csm_contour_map_ce(wks,TEMP2M(i,:,:),res)

  attach_shapefile_outlines(wks,plot)

  draw(plot)

  frame(wks)

  i=i+1

end do

end

also the suprocedure generate two outpus wks and map but i don''t use
nowhere map

can someone help telling me where put the hands for a fast output?

Regards

 ugo
                                               

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Aug 15 11:44:23 2011