map projection problem in pyngl

   [Note to Ufuk Utku Turuncoglu: I think you posted your question
    to pyngl-talk before your registration request had been officially
    okayed, so your posting was automatically rejected. I have since
    okayed your request and am re-posting your question along with some
    answers.]

> hi,
>
> I write a script to plot but it gives following error,
>
> warning:wkColorMap is not a valid resource in contour at this time
> fatal:NhlCvtGenArrayToScalar:GenArray to Scalar with a 0 size array

I believe this is due to the fact that you are reusing a resource list
("res") that you used to set a colormap to then set some contour
resources. The contour plotting function doesn't recognize the "wk"
resources.

What you want to do here is either delete "res" after you set the
colormap, and then create a new "res" resource list, or else rename
the workstation resource list to something else, like "wkres". The
easiest is to probably delete the resource list and create a new one:

res = Ngl.Resources()
#res.wkColorMap = cmap
res.wkColorMap = ["gui_default"]
#res.wkColorMap = ["wgne15"]
#res.wkColorMap = ["rainbow"]
Ngl.set_values(wks, res)
del res

res = Ngl.Resources() # Create new resource list.
res.nglSpreadColors = True
...

> warning:Error retrieving resource mpRightCornerLatF from args - Ignoring Arg
> fatal:NhlCvtGenArrayToScalar:GenArray to Scalar with a 0 size array

I think these error messages have something to do with the way Numeric
handles its "scalar" values. We've had problems in the past, especially
with the upgrade from Numeric 23 to 24, in which scalars need to be
treated differently.

While I'm not sure I understand why this works, try coercing these
values to float:

res.mpLeftCornerLatF = float(lat[0][58])
res.mpLeftCornerLonF = float(lon[0][58])
res.mpRightCornerLatF = float(lat[73][0])
res.mpRightCornerLonF = float(lon[73][0])

I *think* this will have the benefit of forcing them to be the right
type that PyNGL is expecting. Meanwhile, I will try to find out why
PyNGL is not accepting these values as is. Can you do me a favor and
print out the type of one of these values and let me know what it
reports?

>
> script is running before upgrading python 2.2 to 2.4 but now it fails. I
> install all python modules again and run the all example successfully.
> but my scripts does not work. I attach my script into the mail, any
> suggestion will be helpful.
>
> best regards
>
> Ufuk Utku Turuncoglu
> Istanbul Technical University
> Informatics Institute

When you upgraded Python, did you also happen to upgrade Numeric? I think
this could very well be the problem. Also, which version of PyNGL are
you running?

Thanks,

--Mary

> from Numeric import *
> from Scientific.IO.NetCDF import *
> from mx.DateTime import *
> from pyclimate.JDTimeHandler import *
> import Ngl
> import sys
>
> def smth9(x,p,q):
> ni = x.shape[0]
> nj = x.shape[1]
> if (ni < 3 or nj < 3):
> print "smth9: both array dimensions must be at least three."
> sys.exit()
> po4 = p/4.
> qo4 = q/4.
>
> output = zeros([ni,nj]).astype(x.typecode())
> for j in xrange(1,nj-1):
> for i in xrange(1,ni-1):
> jm1 = j-1
> jp1 = j+1
> im1 = i-1
> ip1 = i+1
> term1 = po4*(x[im1,j]+x[i,jm1]+x[ip1,j]+x[i,jp1]-4.*x[i,j])
> term2 = qo4*(x[im1,jp1]+x[im1,jm1]+x[ip1,jm1]+x[ip1,jp1]-4.*x[i,j])
> output[i,j] = float(x[i,j]) + term1 + term2
>
> output[0,:] = x[0,:]
> output[ni-1,:] = x[ni-1,:]
> output[:,0] = x[:,0]
> output[:,nj-1] = x[:,nj-1]
> return output
>
> def ranges2list(strranges):
> rval = []
> inputdata = strranges.split(",")
> for token in inputdata:
> if "-" in token:
> inirange, endrange = token.split("-")
> rval.extend(range(int(inirange), int(endrange)+1))
> else:
> rval.append(int(token))
> return rval
>
> # get command line arguments
> c1 = sys.argv[1] # Start Year
> c2 = sys.argv[2] # End Year
> c3 = sys.argv[3] # Winter or Summer season DJF -- JJA
>
> # set arguments to variable
> yrini = int(c1)
> yrend = int(c2)
>
> # print information
> print "Year Interval ", yrini, yrend
> if c3 == "DJF":
> print "Plotting DJF Average Surface Temperature"
> MONTHS = ranges2list("12,1,2")
> ptitle = "DJF Avg. Surface Temperature (%4d-%4d)" % (yrini, yrend)
> elif c3 == "JJA":
> print "Plotting JJA Average Surface Temperature"
> MONTHS = ranges2list("6,7,8")
> ptitle = "JJA Avg. Surface Temperature (%4d-%4d)" % (yrini, yrend)
> else:
> print "Third argument must be DJF or JJA !!!"
> print "Exit plot ..."
> sys.exit()
>
> # mm5 model output
> nc = NetCDFFile("./mmout_tmp.nc", "r")
> tmean = nc.variables["tmean"]
> t = nc.variables["time"]
> lon = nc.variables["longicrs"]
> lat = nc.variables["latitcrs"]
> xc = nc.variables["x_cross"]
> yc = nc.variables["y_cross"]
> jdth = JDTimeHandler(t.units)
>
> # calculate climatology of model output
> count = 0
> clim = zeros(tmean.shape[1:])
> for irec in range(len(tmean)):
> yy, mm = jdth.getdatefields(t[irec], 2)
> if yrini <= yy <= yrend and mm in MONTHS:
> clim = clim + tmean[irec]
> count = count + 1
> clim = clim/float(count)
>
> # smooth boundries of the domain, 3 grid point average
> im, jm = clim.shape
> clim[:,0] = (clim[:,1]+clim[:,2]+clim[:,3])/3.0
> clim[im-1,:] = (clim[im-4,:]+clim[im-3,:]+clim[im-2,:])/3.0
> # smooth whole domain, light smoothing
> clim = smth9(clim, 0.5, -0.25)
>
> # cru precipitation data
> nc_cru = NetCDFFile("./cru_ts_2_10.1901-2002.tmp.nc")
> cru = nc_cru.variables["tmp"]
> lon_cru = nc_cru.variables["lon"]
> lat_cru = nc_cru.variables["lat"]
>
> # calculate climatology of cru
> clim_cru = zeros(cru.shape[1:])
> yy2 = 1900
> mm2 = 0
> count = 0
> for irec in range(len(cru)):
> if (irec%12 == 0):
> yy2 = yy2 + 1
> mm2=(irec+1)%12
> if (mm2 == 0):
> mm2 = 12
> if yrini <= yy2 <= yrend and mm2 in MONTHS:
> clim_cru = clim_cru + cru[irec]
> count = count + 1
>
> # Plotting using NCL
> # 1 - Create colormap
> cmap = array([[255, 255, 255], [ 0, 0, 0], [255, 255, 255], \
> [244, 255, 244], [217, 255, 217], [163, 255, 163], \
> [106, 255, 106], [ 43, 255, 106], [ 0, 224, 0], \
> [ 0, 134, 0]], Float)/255.0
>
> # 2 - Create NCL environment
> #wks_type = "x11"
> wks_type = "ps"
> #wks_type = "pdf"
> outname = "tmp_%s_%d_%d" % (c3, yrini, yrend)
> wks = Ngl.open_wks(wks_type, outname)
> res = Ngl.Resources()
> #res.wkColorMap = cmap
> res.wkColorMap = ["gui_default"]
> #res.wkColorMap = ["wgne15"]
> #res.wkColorMap = ["rainbow"]
> Ngl.set_values(wks, res)
>
> res.nglSpreadColors = True
>
> res.cnFillOn = True
> res.cnLineLabelsOn = False
> res.cnInfoLabelOn = False
> res.cnLinesOn = False
>
> #res.mpDataBaseVersion = "RANGS_GSHHS"
> res.mpProjection = "LambertConformal"
> res.mpGridAndLimbOn = True
> res.mpAreaMaskingOn = True
> res.mpOutlineOn = True
> res.mpFillOn = False
> res.mpLimitMode = "Corners"
> res.mpLeftCornerLatF = lat[0][58]
> res.mpLeftCornerLonF = lon[0][58]
> res.mpRightCornerLatF = lat[73][0]
> res.mpRightCornerLonF = lon[73][0]
> res.mpLambertParallel1F = 60.
> res.mpLambertParallel2F = 30.
> res.mpLambertMeridianF = 26.
> res.mpPerimOn = True
> res.tfDoNDCOverlay = True
>
> res.lbLabelBarOn = False
> res.pmTickMarkDisplayMode = "Always"
> res.pmLabelBarWidthF = 0.1
> res.tiXAxisString = "Longitude"
> res.tiYAxisString = "Latitude"
> res.tiXAxisFontHeightF = 0.012
> res.tiYAxisFontHeightF = 0.013
> res.tiMainFontHeightF = 0.016
> res.tmXBLabelFontHeightF = 0.012
> res.tmYLLabelFontHeightF = 0.012
> res.tmYRLabelsOn = False
> res.tmXTLabelsOn = False
>
> # create panel plot
> res.nglDraw = False
> res.nglFrame = False
>
> # draw plots
> plot = []
>
> # plot 1: model
> clim_tr = transpose(clim)
> plot.append(Ngl.contour_map(wks, clim_tr, res))
>
> # plot 2: cru
> res.sfMissingValueV = clim_cru[0][0]
> res.mpProjection = "CylindricalEquidistant"
> res.mpLimitMode="LatLon"
> res.mpMinLatF = 18.75
> res.mpMaxLatF = 70.25
> res.mpMinLonF = -15.25
> res.mpMaxLonF = 70.25
> ymin = Ngl.ind(lat_cru[:] == 18.75)[0]
> ymax = Ngl.ind(lat_cru[:] == 70.25)[0]
> xmin = Ngl.ind(lon_cru[:] == -15.25)[0]
> xmax = Ngl.ind(lon_cru[:] == 70.25)[0]
> pp = zeros((ymax-ymin, xmax-xmin), Float)
> pp = clim_cru[ymin:ymax,xmin:xmax]
> plot.append(Ngl.contour_map(wks, pp, res))
>
> # panel
> panelres = Ngl.Resources()
> panelres.nglFrame = False
> panelres.lbLabelBarOn = True
> panelres.lbTitleOn = True
> panelres.lbTitleString = "deg C"
> panelres.lbTitlePosition = "Bottom"
> panelres.lbTitleFontHeightF= 0.012
> panelres.lbTitleDirection = "Across"
> panelres.lbLabelPosition = "Right"
> panelres.nglPanelLabelBarWidthF = 0.42
> panelres.nglPanelLabelBarLabelFontHeightF = 0.012
> panelres.nglPanelLabelBarParallelPosF = 0.27
> panelres.nglPanelLabelBar = True
> panelres.nglPanelFigureStrings = ["MM5", "CRU"]
> panelres.nglPanelFigureStringsJust = "BottomRight"
> Ngl.panel(wks, plot, [1,2], panelres)
>
> # title
> textres = Ngl.Resources()
> textres.txFontHeightF = 0.014
> Ngl.text_ndc(wks, ptitle, 0.5, 0.76, textres)
> Ngl.frame(wks)
>
> Ngl.end()
_______________________________________________
pyngl-talk mailing list
pyngl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/pyngl-talk
Received on Thu Aug 10 2006 - 09:31:45 MDT