Re: (no subject)

From: ugo merlini <ugomerlini_at_nyahnyahspammersnyahnyah>
Date: Tue Aug 16 2011 - 03:45:28 MDT

Hi

I take 50 seconds

Is possibile to create a shapefile map only one time an after use it as background for ploting over an input data for example 2m temperature at time t, t+1 (other loop)? I file for each time.

this way should be more fast

now I'm trying with overlay function so before I'm creating two maps one with the shapefile only the other one the contour only here the code but I'm getting thes warning messages

warning:mpLimitMode is not a valid resource in JAL_contour at this time
warning:mpMinLonF is not a valid resource in JAL_contour at this time
warning:mpMaxLonF is not a valid resource in JAL_contour at this time
warning:mpMinLatF is not a valid resource in JAL_contour at this time
warning:mpMaxLatF is not a valid resource in JAL_contour at this time

code

;----------------------------------------------------------------------
; 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 = 35.0
  res@mpMaxLatF = 48.0
  res@mpMinLonF = 6.5
  res@mpMaxLonF = 19.0
  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(0,:,:) ; 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
  
  plot = gsn_csm_contour(wks,TEMP2M,res)
  overlay(map,plot)
  
;---We're done adding primitives, draw everything and advance frame.
  
  draw(map)
  frame(wks)
end

ugo

> Date: Mon, 15 Aug 2011 13:15:22 -0600
> From: dave.allured@noaa.gov
> Subject: Re: (no subject)
> To: ugomerlini@hotmail.com
> CC: ncl-talk@ucar.edu
>
> Drawing shape files in NCL is slow because of nested loops. On my
> computer I get these speeds to make only one plot:
>
> 3.3 million points X11 output 24 seconds
> 3.3 million points PDF output 36 seconds
>
> What speed are you getting?
>
> --Dave
>
> On 8/15/2011 11:44 AM, ugo merlini wrote:
> > i made my first script which use a shapefile as source but it is
> > very slow to produce the output
                                               

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

This archive was generated by hypermail 2.1.8 : Mon Aug 22 2011 - 08:13:37 MDT