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