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

Hi Ugo,

You might experiment with the resource "gsnAddCyclic" to deal with the gap in contouring over the Meridian:

http://www.ncl.ucar.edu/Document/Graphics/Resources/gsn.shtml

As for limiting the range of lat/lon read from the OpenDAP file, check out the so-called "constraint expressions", in opendap parlance:

http://docs.opendap.org/index.php/UserGuideOPeNDAPMessages#Constraint_Expressions

I don't know about the res/res1 question.

Hope some of that helps...
Rick

On Aug 16, 2011, at 7:45 AM, ugo merlini wrote:

> 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

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