Hello NCL-talk,
I am a new user of ncl, then I am sure there are tons I am missing..
I'd like to plot 2-dimensional data in a vertical height- space (columns)
and (thanks to ncl!!) I succeded in doing it with one file.
Now, I'd like to use contour to plot contouring on top of shading and I
fail to do the right transformation of coordinates to plot both using
a logLinPlotClass, see the script below.
Thank you for your time.
Paolina
*******************************************************************************
;================================================;
; prova-114.ncl
;================================================;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; ================================================;
begin
;=================================================;
; open file and read in data
;=================================================;
ncol=25
nlev=40
;
data1=asciiread ("aveM-pt-bin-114.dat",(/nlev,ncol/),"float")
;
u=data1
u!0="lev"
u!1="col"
u&col=ispan(1,25,1)
;
printVarSummary(u)
; print(u)
;---read data
;
data = asciiread("Mave-masshzm-bin-114.dat",(/5,5/),"float")
z = data(:,:)
printVarSummary(z)
nlat=5
nlon=5
hzmAvg=new( (/nlat,nlon/), float )
hzmAvg = z
;
hzmAvg!0="lat"
hzmAvg!1="lon"
;
printVarSummary(hzmAvg)
;--------------------------------------------------------------------
; reorganize hzmAvg in 1-D array hzmAvg1d
; rearrange HzmAvg1D in ascending order, sort hzmAvg1d also
;--------------------------------------------------------------------
hzmAvg1d =new((/25/),float)
hzmAvg1d = ndtooned(hzmAvg)
hzmAvg1d!0 ="col"
ipAvg = dim_pqsort_n(hzmAvg1d, 2,0) ;no reordering needed
;
printVarSummary(ipAvg)
print(ipAvg)
printVarSummary(hzmAvg1d)
u3d=new((/nlev,ncol/),float)
u3d!0="lev"
u3d!1="col"
u3d=u(:,ipAvg)
u3d&col=u&col(ipAvg)
printVarSummary(u3d)
lev=(/-25.00,25.00,76.98,133.38,195.20,263.67,340.25,426.69,525.07,637.83,767.79,918.14,1092.45,1294.61,1528.71,1798.91,2109.22,2463.31,2864.20,3314.08,3814.08,4364.20,4963.31,5609.22,6298.91,7028.71,7794.61,8592.45,9418.14,10267.79,11137.83,12025.07,12926.69,13840.24,14763.67,15695.20/)
nlev=36
ncol=25
nc=ncol-1
p = new ((/nlev,ncol/), "float" )
p(0,:)=lev(0)
p(1,:)=lev(1)
p(2,:)=lev(2)
p(3,:)=lev(3)
p(4,:)=lev(4)
p(5,:)=lev(5)
p(6,:)=lev(6)
p(7,:)=lev(7)
p(8,:)=lev(8)
p(9,:)=lev(9)
p(10,:)=lev(10)
p(11,:)=lev(11)
p(12,:)=lev(12)
p(13,:)=lev(13)
p(14,:)=lev(14)
p(15,:)=lev(15)
p(16,:)=lev(16)
p(17,:)=lev(17)
p(18,:)=lev(18)
p(19,:)=lev(19)
p(20,:)=lev(20)
p(21,:)=lev(21)
p(22,:)=lev(22)
p(23,:)=lev(23)
p(24,:)=lev(24)
p(25,:)=lev(25)
p(26,:)=lev(26)
p(27,:)=lev(27)
p(28,:)=lev(28)
p(29,:)=lev(29)
p(30,:)=lev(30)
p(31,:)=lev(31)
p(32,:)=lev(32)
p(33,:)=lev(33)
p(34,:)=lev(34)
p(35,:)=lev(35)
; print(p)
printVarSummary(p)
u2 = new ((/36,ncol/), typeof(u3d) )
u2= u3d(0:35,:)
u2!0="lev"
u2!1="col"
du2 =new((/2,36,ncol/),"float")
du2!0="var"
du2!1="lev"
du2!2="col"
du2v = conform_dims(dimsizes(du2), u2,(/1,2/))
du2v(0,:,:) = p
; print(du2v(1,:,:))
printVarSummary(du2v)
;
dims1=(/nlev,ncol/)
range1= new(dims1,"integer")
; printVarSummary(range1)
ip=1
;
do j=0,(nlev-1)
do i= 0, nc
range1(j,i)=ip
;
if(ip.eq.25)then
ip=0
end if
;
ip=ip+1
end do
end do
; print(range1)
;
printVarSummary(range1)
dims = (/2,nlev,ncol/)
coords = new(dims,"integer")
coords(0,:,:) = floattointeger(p)
coords(1,:,:) = range1
coords!0 = "vcoords" ; name dimensions so we can reorder
coords!1 = "lev"
coords!2 = "col"
; print(coords(1,:,:))
printVarSummary(coords)
;
d = du2v(1,:,:)
;********************************
; second plot *
;********************************
data2=asciiread ("aveM-massphow-bin-114.dat",(/nlev,ncol/),"float")
;
v=data2
v!0="lev"
v!1="col"
v&col=ispan(1,25,1)
;
printVarSummary(v)
; print(v)
;*****************************************
; rearrange data according to ipAvg
;*****************************************
v3d=new((/nlev,ncol/),float)
v3d!0="lev"
v3d!1="col"
v3d=v(:,ipAvg)
v3d&col=v&col(ipAvg)
printVarSummary(v3d)
v2 = new ((/36,ncol/), typeof(v3d) )
v2= v3d(0:35,:)
v2!0="lev"
v2!1="col"
dv2 =new((/2,36,ncol/),"float")
dv2!0="var"
dv2!1="lev"
dv2!2="col"
dv2v = conform_dims(dimsizes(dv2), v2,(/1,2/))
dv2v(0,:,:) = p
; print(dv2v(1,:,:))
printVarSummary(dv2v)
;
dims1=(/nlev,ncol/)
range2= new(dims1,"integer")
; printVarSummary(range2)
ip=1
;
do j=0,(nlev-1)
do i= 0, nc
range2(j,i)=ip
;
if(ip.eq.25)then
ip=0
end if
;
ip=ip+1
end do
end do
; print(range2)
;
printVarSummary(range2)
dims = (/2,nlev,ncol/)
coords2 = new(dims,"integer")
coords2(0,:,:) = floattointeger(p)
coords2(1,:,:) = range2
coords2!0 = "vcoords" ; name dimensions so we can reorder
coords2!1 = "lev"
coords2!2 = "col"
; print(coords2(1,:,:))
printVarSummary(coords2)
dd = dv2v(1,:,:)
;**************************************************
; read in coordinates
;**************************************************
lev1 = coords(vcoords|0,lev|:,col|:)
col = coords(vcoords|1,lev|:,col|:)
; the challenge with this type of data is a two-dimensional vertical
; coordinate. The gsn_csm high-level graphical interfaces can handle 2D
; lat/lon coordinates but are not currently adpated for 2D vertical coords.
lev2 = coords2(vcoords|0,lev|:,col|:)
col2 = coords2(vcoords|1,lev|:,col|:)
;**************************************************
; assign variable meta data
;**************************************************
wks = gsn_open_wks("ps","psi-pt-114") ; open a ps file
gsn_define_colormap(wks,"gui_default") ; choose colormap
i = NhlNewColor(wks,0.8,0.8,0.8) ; add a gray to colormap
res = True ; plot mods desired
; res@gsnDraw = False ; don't draw
res@gsnFrame = False ; don't advance frame
res@sfXArray = col ; could be reduced to 1D
res@sfYArray = lev1 ; 2D
res@cnInfoLabelOn = False ; turn off contour info label
res@cnLevelSelectionMode = "ExplicitLevels" ; manual levels
res@cnMinLevelValF = 0.
res@cnMaxLevelValF = 400.
res@cnLevelSpacingF = 5.
res@cnFillOn = True ; turn on color
res@cnLineLabelsOn = False ; turn off line labels
res@gsnSpreadColors = True ; use full range of colormap
res@cnLinesOn = False ; turn off contour lines
; res@cnLabelMasking = True
; Linearize the plot by overlaying on a logLinPlot. It's important to use the
; "curvilinear" gridType instead of the default 2D "spherical" grid type
; because the spherical grid type assumes that the X coordinates are modular
; and that the Y coordinates only range from -90 to 90.
plot = gsn_contour(wks,d,res) ; contour the variable
;**************************************************
;
res1 = res
res1@cnFillOn = False ; turn on color
res1@cnLinesOn = True ; turn off contour lines
res@sfXArray = col ; could be reduced to 1D
res@sfYArray = lev1 ; 2D
; delete(res1@cnLevels)
; res1@cnLevelSelectionMode = "EqualSpacedLevels" ; set equal spaced
contour levels
; res1@cnMaxLevelCount = 20
res1@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels
res1@cnLevels =
(/-4.E-03,-3.5E-03,-3.E-03,-2.5E-03,-2.E-03,-1.5E-03,-1.E-03,-.5E-02,1.E-03,
\
1.5E-03,2.E-03,2.5E-03,3.E-03,3.5E-03,4.E-03,4.5E-03,5.E-03,5.5E-03,6.E-03/)
; res1@cnMinLevelValF = -4.E-03
; res1@cnMaxLevelValF = 8.E-03
; res1@cnLevelSpacingF = 0.005
res1@gsnContourNegLineDashPattern = 1 ; sets negative contours
to dash pattern 1
res1@gsnContourPosLineDashPattern = 0
; res1@gsnContourZeroLineThicknessF = 2. ; doubles thickness of
zero contour
plotv = gsn_contour(wks,dd,res1) ; contour the variable
;
setvalues plot
"trGridType" : "curvilinear"
"tiMainString" : "str-rrc day 74/76"
"tiXAxisString" : "Column " ; x-axis title
"tiYAxisString" : "Height (m)" ; y-axis title
end setvalues
ll = create "ll" logLinPlotClass wks
"trXMinF" : min(col)
"trXMaxF" : max(col)
"trYMinF" : min(lev1)
"trYMaxF" : max(lev1)
"pmTickMarkDisplayMode" : "always"
end create
;
overlay(ll, plot)
maximize_output(wks,False)
draw ((/plot,plotv/))
frame(wks)
end
***************************************************************************
--
Dr. Paolina Bongioannini Cerlini
University of Perugia
Dip.to Ingegneria
CIRIAF -room 8
Via G. Duranti, 67
06125 - Perugia
Tel.+39 075 585 3576
Fax +39 075 585 3697
e-mail: paolina.cerlini@unipg.it
paolina.cerlini@fisica.unipg.it
cell.:+39 3405761980
--
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 14 09:24:24 2013
This archive was generated by hypermail 2.1.8 : Tue Mar 19 2013 - 16:31:18 MDT