KevinMidboe/linguist
1
0
Fork 0
mirror of https://github.com/KevinMidboe/linguist.git synced 2025-12-08 04:18:48 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into more-encompassing-number-skips

Browse Source
This commit is contained in:
Arfon Smith
2015-07-29 13:54:51 +01:00
parent 885b5aab41
commit 90a293727d
434 changed files with 79876 additions and 56121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

109
samples/NCL/PrnOscPat_driver.ncl Normal file
View File

@@ -0,0 +1,109 @@
undef("PrnOscPat_driver")
function PrnOscPat_driver(eof[*][*][*]:numeric, eof_ts[*][*]:numeric, kPOP[1]:integer)
; =================================================================
; compute Principal Oscillation Patterns (POPs)
; =================================================================
local dim_ts, dim_eof, neof, ntim, nlat, mlon, dnam_ts, dnam_eof, neof, j \
, cov0, cov1, cov0_inverse, A, z, Z, pr, pi, zr, zi, mean, stdev \
, evlr, eigi, eigr
begin
dim_ts = dimsizes(eof_ts) ; (neof,ntim)
dim_eof = dimsizes(eof) ; (neof,nlat,mlon)
ntim = dim_ts(1)
neof = dim_eof(0)
nlat = dim_eof(1)
mlon = dim_eof(2)
dnam_ts = getvardims(eof_ts) ; dimension names
dnam_eof= getvardims(eof) ; used at end for meta data
; =================================================================
; lag-0 and lag-1 matrices
; =================================================================
if (get_ncl_version().eq."6.1.2") then ; bug in 6.1.2
cov0 = covcorm(eof_ts,(/1,0/)) ; lag-0 covariance matrix
else
cov0 = covcorm(eof_ts,(/0,1/)) ; lag-0 covariance matrix (n x n)
end if
; either
cov1 = covcorm_xy(eof_ts, eof_ts, (/0,1,0/)) ; lag-1
;cov1 = covcorm_xy(eof_ts(:,0:ntim-2) \ ; alternative, brute force
; ,eof_ts(:,1:ntim-1), (/0,0,0/))
;printVarSummary(cov1)
; =================================================================
; matrix A contains information for evolution of the POP system.
; POPs are eigenvectors of A.
; =================================================================
cov0_inverse = inverse_matrix(cov0)
A = cov1#inverse_matrix(cov0) ; [*][*] => neof x neof
; =================================================================
; NCL 6.1.1 of dgeevx: evlr(2,2,N,N) ; (left(0)/right(1), real(0)/imag(1),:,:)
; Eigenvalues are returned as attributes: eigi = evlr@eigi ; eigr = evlr@eigr
; =================================================================
evlr = dgeevx_lapack(A, "B", "V", "V", "B", False)
; =================================================================
; POP time series from eigenvalues and right eigenvectors
; =================================================================
;PR = (/ evlr(1,0,:,:) /) ; right ev (1), real part (0)
;PI = (/ evlr(1,1,:,:) /) ; right ev (1), imag part (1)
; kPOP is what we want; use righteigenvector
pr = (/ evlr(1,0,kPOP-1,:) /) ; right ev (1), real part (0), row 'kPOP-1'
pi = (/ evlr(1,1,kPOP-1,:) /) ; right ev (1), imag part (1), row 'kPOP-1'
z = inverse_matrix( (/ (/sum(pr*pr), sum(pr*pi)/) \
, (/sum(pr*pi), sum(pi*pi)/) /))#(/pr,pi/)#eof_ts
; complex conjugate
z = (/z(0,:), -z(1,:)/) ; real & imag series
z = dim_rmvmean_n(z,1)
mean = dim_avg_n(z,1) ; calculate mean
stdev= dim_stddev_n(z,1) ; calculate stdev
z = dim_standardize_n(z,1,1) ; standardize time series
z!0 = "nPOP" ; add meta data
z!1 = dnam_ts(1)
z&nPOP = (/0,1/)
z&$dnam_ts(1)$ = eof_ts&$dnam_ts(1)$
z@stdev = stdev
z@mean = mean
z@long_name = "POP timeseries"
;printVarSummary(z)
; =================================================================
; POP spatial patterns
; =================================================================
zr = pr(0)*eof(0,:,:) ; construct POP spatial domain
zi = pi(0)*eof(0,:,:)
do j=1,neof-1
zr = zr + pr(j)*eof(j,:,:)
zi = zi + pi(j)*eof(j,:,:)
end do
Z = (/zr*stdev(0), -zi*stdev(1)/) ; scale patterns by time series stdev
Z!0 = "nPOP" ; add meta data
Z!1 = dnam_eof(1)
Z!2 = dnam_eof(2)
Z&nPOP = (/0,1/)
Z&$dnam_eof(1)$ = eof&$dnam_eof(1)$
Z&$dnam_eof(2)$ = eof&$dnam_eof(2)$
Z@long_name = "POP pattern"
;printVarSummary(Z)
; =================================================================
; return POP time series and POP spatial patterns as a
; variable of type 'list' which contains 2 variables
; =================================================================
return( [/z, Z/] ) ; this is type "list"
end

115
samples/NCL/WRF_static_2.ncl Normal file
View File

@@ -0,0 +1,115 @@
;*************************************************
; WRF static: panel different variables
;************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
begin
;************************************************
; open file and read in data
;************************************************
f = addfile("static.wrfsi.nc", "r")
;************************************************
; Read variables
;************************************************
use = f->use(0,0,:,:) ; land use dominant category
stl = f->stl(0,0,:,:) ; top layer (0-30cm) dom cat soiltype
sbl = f->sbl(0,0,:,:) ; bottom layer (30-90cm) dom cat soiltype
lat2d = f->lat(0,0,:,:)
lon2d = f->lon(0,0,:,:)
lsMask= f->lnd(0,0,:,:) ; land (1) water (0) mas
;************************************************
; Use mask function to set all ocean areas to _FillValue
;************************************************
use = mask(use,lsMask,1)
stl = mask(stl,lsMask,1)
sbl = mask(sbl,lsMask,1)
;************************************************
; Associate 2D coordinates with variables for plotting
;************************************************
use@lat2d = lat2d
use@lon2d = lon2d
stl@lat2d = lat2d
stl@lon2d = lon2d
sbl@lat2d = lat2d
sbl@lon2d = lon2d
;************************************************
; The file should be examined via: ncdump -v grid_type static.wrsi
; This will print the print type. then enter below.
;************************************************
projection = "mercator"
;************************************************
; create plots
;************************************************
wks = gsn_open_wks("ps" ,"WRF_static") ; ps,pdf,x11,ncgm,eps
gsn_define_colormap(wks ,"BlAqGrYeOrReVi200"); choose colormap
res = True ; plot mods desired
res@gsnSpreadColors = True ; use full range of colormap
res@cnFillOn = True ; color plot desired
res@cnLinesOn = False ; turn off contour lines
res@cnLineLabelsOn = False ; turn off contour labels
res@cnLevelSpacingF = 1 ; manually specify interval
res@cnFillMode = "RasterFill" ; activate raster mode
res@lbLabelAutoStride = True ; let NCL figure lb stride
;************************************************
; Turn on lat / lon labeling
;************************************************
;;res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
res@mpProjection = projection
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat2d(0,0)
res@mpLeftCornerLonF = lon2d(0,0)
res@mpRightCornerLatF = lat2d(nlat-1,mlon-1)
res@mpRightCornerLonF = lon2d(nlat-1,mlon-1)
res@mpCenterLonF = f->LoV ; set center logitude
if (projection.eq."LambertConformal") then
res@mpLambertParallel1F = f->Latin1
res@mpLambertParallel2F = f->Latin2
res@mpLambertMeridianF = f->LoV
end if
res@mpFillOn = False ; turn off map fill
res@mpOutlineDrawOrder = "PostDraw" ; draw continental outline last
res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
;;res@tfDoNDCOverlay = True ; True only for 'native' grid
res@gsnAddCyclic = False ; data are not cyclic
;************************************************
; allocate array for 3 plots
;************************************************
plts = new (3,"graphic")
;************************************************
; Tell NCL not to draw or advance frame for individual plots
;************************************************
res@gsnDraw = False ; (a) do not draw
res@gsnFrame = False ; (b) do not advance 'frame'
plts(0) = gsn_csm_contour_map(wks,use,res)
plts(1) = gsn_csm_contour_map(wks,stl,res)
plts(2) = gsn_csm_contour_map(wks,sbl,res)
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
resP = True ; modify the panel plot
resP@txString = "Land Use and Soil Type"
resP@gsnMaximize = True ; maximize panel area
resP@gsnPanelRowSpec = True ; specify 1 top, 2 lower level
gsn_panel(wks,plts,(/1,2/),resP) ; now draw as one plot
end

160
samples/NCL/WRF_track_1.ncl Normal file
View File

@@ -0,0 +1,160 @@
;********************************************************
; Plot storm stracks from wrfout files.
;********************************************************
;
; JUN-18-2005
; So-Young Ha (MMM/NCAR)
; SEP-01-2006
; Slightly modified by Mary Haley to add some extra comments.
; ===========================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
; DATES
date = (/1512,1600,1612,1700,1712,1800,1812,1900/)
ndate = dimsizes(date)
sdate = sprinti("%4.0i",date)
; Experiment name (for legend)
EXP = (/"EXP_I"/) ; (/"EXP_I","EXP_II","EXP_III"/)
nexp = dimsizes(EXP)
; To get lat/lon info.
a = addfile("wrfout_d01_2003-07-15_00:00:00.nc","r")
lat2d = a->XLAT(0,:,:)
lon2d = a->XLONG(0,:,:)
dimll = dimsizes(lat2d)
nlat = dimll(0)
mlon = dimll(1)
; Sea Level Pressure
slp = wrf_user_getvar(a,"slp",0)
dims = dimsizes(slp)
; Array for track
time = new(ndate,string)
imin = new(ndate,integer)
jmin = new(ndate,integer)
smin = new(ndate,integer)
; =======
; ndate
; =======
fs = systemfunc("ls wrfout*00")
nfs= dimsizes(fs)
if(nfs .ne. ndate) then
print("Check input data:"+nfs+" .ne. "+ndate)
end if
do ifs=0,nfs-1
f = addfile(fs(ifs)+".nc","r")
time(ifs) = wrf_user_list_times(f)
; print(time(ifs))
slp2d = wrf_user_getvar(f,"slp",0)
; We need to convert 2-D array to 1-D array to find the minima.
slp1d = ndtooned(slp2d)
smin(ifs) = minind(slp1d)
; Convert the index for 1-D array back to the indeces for 2-D array.
minij = ind_resolve(ind(slp1d.eq.min(slp2d)),dims)
imin(ifs) = minij(0,0)
jmin(ifs) = minij(0,1)
; print(time(ifs)+" : "+min(slp2d)+" ("+imin(ifs)+","+jmin(ifs)+")")
end do
;
; Graphics section
wks=gsn_open_wks("ps","track") ; Open PS file.
gsn_define_colormap(wks,"BlGrYeOrReVi200") ; Change color map.
res = True
res@gsnDraw = False ; Turn off draw.
res@gsnFrame = False ; Turn off frame advance.
res@gsnMaximize = True ; Maximize plot in frame.
res@tiMainString = "Hurricane Isabel" ; Main title
WRF_map_c(a,res,0) ; Set up map resources
; (plot options)
plot = gsn_csm_map(wks,res) ; Create a map.
; Set up resources for polymarkers.
gsres = True
gsres@gsMarkerIndex = 16 ; filled dot
;gsres@gsMarkerSizeF = 0.005 ; default - 0.007
cols = (/5,160,40/)
; Set up resources for polylines.
res_lines = True
res_lines@gsLineThicknessF = 3. ; 3x as thick
dot = new(ndate,graphic) ; Make sure each gsn_add_polyxxx call
line = new(ndate,graphic) ; is assigned to a unique variable.
; Loop through each date and add polylines to the plot.
do i = 0,ndate-2
res_lines@gsLineColor = cols(0)
xx=(/lon2d(imin(i),jmin(i)),lon2d(imin(i+1),jmin(i+1))/)
yy=(/lat2d(imin(i),jmin(i)),lat2d(imin(i+1),jmin(i+1))/)
line(i) = gsn_add_polyline(wks,plot,xx,yy,res_lines)
end do
lon1d = ndtooned(lon2d)
lat1d = ndtooned(lat2d)
; Loop through each date and add polymarkers to the plot.
do i = 0,ndate-1
print("dot:"+lon1d(smin(i))+","+lat1d(smin(i)))
gsres@gsMarkerColor = cols(0)
dot(i)=gsn_add_polymarker(wks,plot,lon1d(smin(i)),lat1d(smin(i)),gsres)
end do
; Date (Legend)
txres = True
txres@txFontHeightF = 0.015
txres@txFontColor = cols(0)
txid1 = new(ndate,graphic)
; Loop through each date and draw a text string on the plot.
do i = 0, ndate-1
txres@txJust = "CenterRight"
ix = smin(i) - 4
print("Eye:"+ix)
if(i.eq.1) then
txres@txJust = "CenterLeft"
ix = ix + 8
end if
txid1(i) = gsn_add_text(wks,plot,sdate(i),lon1d(ix),lat1d(ix),txres)
end do
; Add marker and text for legend. (Or you can just use "pmLegend" instead.)
txres@txJust = "CenterLeft"
txid2 = new(nexp,graphic)
pmid2 = new(nexp,graphic)
do i = 0,nexp-1
gsres@gsMarkerColor = cols(i)
txres@txFontColor = cols(i)
ii = ((/129,119,109/)) ; ilat
jj = ((/110,110,110/)) ; jlon
ji = ii*mlon+jj ; col x row
pmid2(i) = gsn_add_polymarker(wks,plot,lon1d(ji(i)),lat1d(ji(i)),gsres)
txid2(i) = gsn_add_text(wks,plot,EXP(i),lon1d(ji(i)+5),lat1d(ji(i)),txres)
end do
draw(plot)
frame(wks)
end

129
samples/NCL/cru_8.ncl Normal file
View File

@@ -0,0 +1,129 @@
;*****************************************************
; cru_8.ncl
;
; Concepts illustrated:
; - Plotting CRU (Climate Research Unit)/ BADC data
; - Selecting a sub-period
; - calculating a climatology
; - Drawing raster contours; very basic graphics
;
;*****************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" ; not needed 6.20 onward
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
; create references (pointers) to the files
diri = "./"
fcld = addfile(diri+"cru_ts3.21.1901.2012.cld.dat.nc", "r")
fdtr = addfile(diri+"cru_ts3.21.1901.2012.dtr.dat.nc", "r")
ffrs = addfile(diri+"cru_ts3.21.1901.2012.frs.dat.nc", "r")
fpet = addfile(diri+"cru_ts3.21.1901.2012.pet.dat.nc", "r")
fpre = addfile(diri+"cru_ts3.21.1901.2012.pre.dat.nc", "r")
ftmn = addfile(diri+"cru_ts3.21.1901.2012.tmn.dat.nc", "r")
ftmp = addfile(diri+"cru_ts3.21.1901.2012.tmp.dat.nc", "r")
ftmx = addfile(diri+"cru_ts3.21.1901.2012.tmx.dat.nc", "r")
fvap = addfile(diri+"cru_ts3.21.1901.2012.vap.dat.nc", "r")
fwet = addfile(diri+"cru_ts3.21.1901.2012.wet.dat.nc", "r")
; specify start & last dates (arbitrary)
ymStrt = 199101
ymLast = 200012
; get index values of start/lat dates
time = fcld->time
yyyymm = cd_calendar(time, -1)
ntStrt = ind(yyyymm.eq.ymStrt) ; index values
ntLast = ind(yyyymm.eq.ymLast)
; read time segment
cld = fcld->cld(ntStrt:ntLast,:,:)
dtr = fdtr->dtr(ntStrt:ntLast,:,:)
frs = ffrs->frs(ntStrt:ntLast,:,:)
pet = fpet->pet(ntStrt:ntLast,:,:)
pre = fpre->pre(ntStrt:ntLast,:,:)
tmn = ftmn->tmn(ntStrt:ntLast,:,:)
tmp = ftmp->tmp(ntStrt:ntLast,:,:)
tmx = ftmx->tmx(ntStrt:ntLast,:,:)
vap = fvap->vap(ntStrt:ntLast,:,:)
wet = fwet->wet(ntStrt:ntLast,:,:)
printVarSummary(cld) ; [time | 120] x [lat | 360] x [lon | 720]
; calculate monthly climatologies
cldclm = clmMonTLL(cld)
dtrclm = clmMonTLL(dtr)
frsclm = clmMonTLL(frs)
petclm = clmMonTLL(pet)
preclm = clmMonTLL(pre)
tmnclm = clmMonTLL(tmn)
tmpclm = clmMonTLL(tmp)
tmxclm = clmMonTLL(tmx)
vapclm = clmMonTLL(vap)
wetclm = clmMonTLL(wet)
printVarSummary(cldclm) ; [month | 12] x [lat | 360] x [lon | 720]
;************************************
; create plots ... very simple
;************************************
nt = 6
month = "July"
yrStrt = ymStrt/100
yrLast = ymLast/100
title = month+": "+yrStrt+"-"+yrLast
wks = gsn_open_wks("ps","cru") ; open a ps file
gsn_define_colormap(wks,"ncl_default") ; choose colormap; not needed 6.20 onward
plot = new(2,graphic) ; create graphic array
res = True
res@cnFillOn = True ; turn on color fill; not needed 6.20 onward
res@cnFillMode = "RasterFill" ; Raster Mode
res@cnLinesOn = False ; Turn off contour lines
res@gsnDraw = False ; do not draw picture
res@gsnFrame = False ; do not advance frame
res@lbOrientation = "Vertical" ; vertical label bar
resp = True
resp@gsnMaximize = True ; make ps, eps, pdf large
resp@txString = title+": CLD, FRS"
plot(0)=gsn_csm_contour_map_ce(wks,cldclm(nt,:,:),res)
plot(1)=gsn_csm_contour_map_ce(wks,frsclm(nt,:,:),res)
gsn_panel(wks,plot,(/2,1/),resp)
resp@txString = title+": PET, VAP"
plot(0)=gsn_csm_contour_map_ce(wks,petclm(nt,:,:),res)
plot(1)=gsn_csm_contour_map_ce(wks,vapclm(nt,:,:),res)
gsn_panel(wks,plot,(/2,1/),resp)
resp@txString = title+": TMN, TMX"
plot(0)=gsn_csm_contour_map_ce(wks,tmnclm(nt,:,:),res)
plot(1)=gsn_csm_contour_map_ce(wks,tmxclm(nt,:,:),res)
gsn_panel(wks,plot,(/2,1/),resp)
resp@txString = title+": TMP, DTR"
plot(0)=gsn_csm_contour_map_ce(wks,tmpclm(nt,:,:),res)
plot(1)=gsn_csm_contour_map_ce(wks,dtrclm(nt,:,:),res)
gsn_panel(wks,plot,(/2,1/),resp)
resp@txString = title+": WET, PRE"
plot(0)=gsn_csm_contour_map_ce(wks,wetclm(nt,:,:),res)
;colors = (/ ... /)
;res@cnFillPalette = colors ; optional: distinct colors for categories
res@cnLevelSelectionMode = "ExplicitLevels" ; use unequal spacing
res@cnLevels = (/2.0,10,25,37.5,50,75,100,125,150,175,200,300,400,500,750/)
plot(1)=gsn_csm_contour_map_ce(wks,preclm(nt,:,:),res)
gsn_panel(wks,plot,(/2,1/),resp)

20
samples/NCL/gsn_csm_xy2_time_series_inputs.ncl Normal file
View File

@@ -0,0 +1,20 @@
;******************** Inputs Regarding Input and Output Data *************************************
;netCDFFilePath = "NULL-MYD04_L2.051-MIL2ASAE.0022-AERONET_AOD_L2.2-20112106165049.nc"
;outputFilePath = "plot-output"
;******************* Inputs Regarding Data Structure ***********************************************
;lPlotVariablesList = "mean_AERONET_AOD_L2_2_AOD0558intrp_Ames,mean_MIL2ASAE_0022_AOD0866b_Ames"
;rPlotVariablesList = "medn_MYD04_L2_051_AOD0550dpbl_l_Ames"
;xDimName = "time"
;xDimSize = 365
;******************* Inputs Regarding the View Annotations ****************************************
;title = "MAPSS Time Series"
;yLAxisLabel = "Mean AOD"
;yRAxisLabel = "Median AOD"
;*******************END INPUTS ********************************************************************

128
samples/NCL/hdf4sds_7.ncl Normal file
View File

@@ -0,0 +1,128 @@
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;**************************************************************
; User Input
;***************************************************************
; INPUT
diri = "./" ; input directory
fili = "wv_LV3_MET08_20050102_12345678_L00013712E00013712.hdf"
pltDir = "./" ; directory for plot output
sfx = get_file_suffix(fili,1)
;pltName = sfx@fBase ; output graphic name
pltName = "hdf4sds"
pltType = "ps"
;***************************************************************
; End User Input
;***************************************************************
;***************************************************************
; Open SEVIRI L3 'wv' HDF file
;***************************************************************
; Note the rather unusual data format: flag *prepended* to data value
;***************************************************************
; integer twc_lv3 ( fakeDim0, fakeDim1 )
; long_name : total water vapour column + flag
; units : fmmmm
; format : I4
; valid_range : ( 10000, 38000 )
; _FillValue : -99
; legend_01 : f = flag
; legend_02 : f = 1 averaged level 2 values
; legend_03 : f = 2 interpolated from averaged level 2 values
; legend_04 : f = 3 gaps filled with NVAP climatology
; legend_05 : mmmm = water vapour column in mm * 100. as integer
; legend_06 : Example: 11025 means: flag = 1, 10.25 mm water vapour column
; min_lat : -74.75
; max_lat : 61.75
; min_lon : -75.25
; max_lon : 75.25
; dlat : 0.5
; dlon : 0.5
;---------------------------------------------------------------
f = addfile (diri+fili, "r")
ifx = f->twc_lv3 ; fmmmm (integer)
printVarSummary(ifx)
flag = ifx/10000 ; extract flag
ix = ifx - flag*10000 ; extract mmmm
x = ix*0.01 ; scale
; create meta data for 'x'
dimx = dimsizes(x)
nlat = dimx(0) ; grid size x(nlat,mlon)
mlon = dimx(1)
lat = fspan(ifx@min_lat, ifx@max_lat, nlat)
lat@units = "degrees_north"
lon = fspan(ifx@min_lon, ifx@max_lon, mlon)
lon@units = "degrees_east"
x!0 = "lat"
x!1 = "lon"
x&lat = lat
x&lon = lon
x@long_name = "SEVIRI: Total Water Vapor"
x@units = "mm"
delete( [/ifx, ix/] ) ; no longer needed
;***************************************************************
; Create plot
;***************************************************************
wks = gsn_open_wks(pltType, pltDir+pltName)
plot = new (2, "graphic")
res = True ; plot mods desired
res@gsnAddCyclic = False ; data noty global
res@gsnDraw = False
res@gsnFrame = False
res@cnFillOn = True ; turn on color fill
res@cnLinesOn = False ; turn of contour lines
res@cnFillMode = "RasterFill" ; Raster Mode
res@cnLinesOn = False ; Turn off contour lines
res@cnLineLabelsOn = False ; Turn off contour lines
res@cnMissingValFillColor= "background" ; "foreground"
res@mpCenterLonF = 0.5*(min(x&lon) + max(x&lon))
res@mpMinLatF = min(x&lat)
res@mpMaxLatF = max(x&lat)
res@mpMinLonF = min(x&lon)
res@mpMaxLonF = max(x&lon)
;res@lbOrientation = "Vertical"
plot(0) = gsn_csm_contour_map_ce(wks,x, res)
; plot flag
copy_VarCoords(x, flag)
flag@long_name = "Flag"
flag@units = "1=avg(L2), 2=int(L2), 3=NVAP"
print(flag&lat+" "+flag(:,{30}))
res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
res@cnMinLevelValF = 2 ; set min contour level
res@cnMaxLevelValF = 3 ; one less than max
res@cnLevelSpacingF = 1 ; set contour spacing
res@lbLabelStrings = ispan(1,3,1) ; 1, 2, 3
res@lbLabelPosition = "Center" ; label position
res@lbLabelAlignment = "BoxCenters"
res@gsnLeftString = ""
res@gsnRightString = ""
res@gsnCenterString = "flag: 1=avg(L2), 2=int(L2), 3=NVAP"
plot(1) = gsn_csm_contour_map_ce(wks,flag, res)
resP = True ; modify the panel plot
resP@txString = fili
resP@gsnMaximize = True
gsn_panel(wks,plot,(/1,2/),resP) ; now draw as one plot

125
samples/NCL/mask_12.ncl Normal file
View File

@@ -0,0 +1,125 @@
;----------------------------------------------------------------------
; mask_12.ncl
;
; Concepts illustrated:
; - Using a worldwide shapefile to create a land/ocean mask
; - Masking a data array based on a geographical area
; - Attaching shapefile polylines to a map plot
; - Attaching lat/lon points to a map using gsn_coordinates
;----------------------------------------------------------------------
; Downloaded GSHHS shapefiles from:
;
; http://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/
;
; Used the "coarsest" one: "GSHHS_shp/c/GSHHS_c_L1.shp".
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "./shapefile_mask_data.ncl"
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
WRITE_MASK = True
DEBUG = False
;---Read data to plot and mask
dir = "$NCARG_ROOT/lib/ncarg/data/cdf/"
cdf_prefix = "uv300"
cdf_file = dir + cdf_prefix + ".nc"
fin = addfile(cdf_file,"r")
u = fin->U(1,:,:)
;
; Create a mask array the same size as "u", using
; lat/lon data read off a shapefile.
;
shpfile = "GSHHS_shp/c/GSHHS_c_L1.shp"
opt = True
opt@return_mask = True
land_mask = shapefile_mask_data(u,shpfile,opt)
;---Mask "u" against land and ocean.
u_land_mask = where(land_mask.eq.1,u,u@_FillValue)
u_ocean_mask = where(land_mask.eq.0,u,u@_FillValue)
copy_VarMeta(u,u_land_mask)
copy_VarMeta(u,u_ocean_mask)
;---Start the graphics
wks = gsn_open_wks("ps","mask")
res = True
res@gsnMaximize = True ; maximize plot in frame
res@gsnDraw = False ; don't draw plot yet
res@gsnFrame = False ; don't advance frame yet
res@cnFillOn = True
res@cnLineLabelsOn = False
res@cnLinesOn = False
;---Make sure both plots have same contour levels
mnmxint = nice_mnmxintvl(min(u),max(u),25,False)
res@cnLevelSelectionMode = "ManualLevels"
res@cnMinLevelValF = mnmxint(0)
res@cnMaxLevelValF = mnmxint(1)
res@cnLevelSpacingF = mnmxint(2)
res@lbLabelBarOn = False
res@gsnAddCyclic = False
res@mpFillOn = False
res@mpOutlineOn = False
res@gsnRightString = ""
res@gsnLeftString = ""
;---Create plot of original data and attach shapefile outlines
res@tiMainString = "Original data with shapefile outlines"
map_data = gsn_csm_contour_map(wks,u,res)
dum1 = gsn_add_shapefile_polylines(wks,map_data,shpfile,False)
;---Create plots of masked data
res@tiMainString = "Original data masked against land"
map_land_mask = gsn_csm_contour_map(wks,u_land_mask,res)
res@tiMainString = "Original data masked against ocean"
map_ocean_mask = gsn_csm_contour_map(wks,u_ocean_mask,res)
if(DEBUG) then
mkres = True
; mkres@gsMarkerSizeF = 0.007
mkres@gsnCoordsAttach = True
gsn_coordinates(wks,map_data,u,mkres)
mkres@gsnCoordsNonMissingColor = "yellow"
mkres@gsnCoordsMissingColor = "black"
gsn_coordinates(wks,map_land_mask,u_land_mask,mkres)
gsn_coordinates(wks,map_ocean_mask,u_ocean_mask,mkres)
end if
;---Add shapefile outlines
dum2 = gsn_add_shapefile_polylines(wks,map_land_mask,shpfile,False)
dum3 = gsn_add_shapefile_polylines(wks,map_ocean_mask,shpfile,False)
;---Draw all three plots on one page
pres = True
pres@gsnMaximize = True
pres@gsnPanelLabelBar = True
gsn_panel(wks,(/map_data,map_land_mask,map_ocean_mask/),(/3,1/),pres)
if(WRITE_MASK) then
delete(fin) ; Close file before we open again.
;
; Make copy of file so we don't overwrite original.
; This is not necessary, but it's safer.
;
new_cdf_file = cdf_prefix + "_with_mask.nc"
system("/bin/cp " + cdf_file + " " + new_cdf_file)
finout = addfile(new_cdf_file,"w")
filevardef(finout, "land_mask", typeof(land_mask), (/ "lat", "lon" /) )
finout->land_mask = (/land_mask/)
end if
end

115
samples/NCL/mcsst_1.ncl Normal file
View File

@@ -0,0 +1,115 @@
;*****************************************************
; mcsst_1.ncl
;
; Concepts illustrated:
; - Plotting NAVO MCSST data
; - Using fbindirread to read in fortran binary data
; - Converting "byte" data to "float"
; - Adding meta data (attributes and coordinates) to a variable
; - Adding gray to an existing color map
; - Spanning all but the last two colors in a color map for contour fill
; - Drawing raster contours
;
;*****************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;***************************************
; type of data available on file
;***************************************
; ipar=0 Weekly Binned Sea Surface Temperature
; ipar=1 Number of Points in Bin
; ipar=2 Weekly Binned Sea Surface Temperature Anomaly
; ipar=3 Interpolated Sea Surface Temperature
; ipar=4 Interpolated Sea Surface Temperature Anomaly
;***************************************
begin
ipar = 3
fname = "2001311d18N16.dat"
tmp = fbindirread(fname,ipar,(/1024,2048/),"byte")
;***************************************
; convert to float and then change to true SST
;***************************************
xslope = 0.15
if(ipar.eq.4.or.ipar.eq.2)then ; anom has different intercept
yint = -20.0
end if
if(ipar.eq.3.or.ipar.eq.0)then
yint = -3.0
end if
sst = new((/1024,2048/),"float") ; create float var
sst = tmp*xslope+yint ; convert to float
delete(tmp) ; delete unecessary array
;***************************************
; assign missing values. The original missing value was zero, but since it was
; not assigned in NCL, it was not recognized. The new missing values are
; listed below. These will be changed later.
;***************************************
if(ipar.eq.4)then
sst@_FillValue = -20
end if
if(ipar.eq.3.or.ipar.eq.0)then
sst@_FillValue = -3
end if
;***************************************
; create coordinate variables
;***************************************
nlat = 1024
dy = 180./nlat
lat = (90. -(ispan(0,1023,1)*dy))-dy/2
lat!0 = "lat"
lat&lat = lat
lat@units = "degrees_north"
nlon = 2048
dx = 360./nlon
lon = (ispan(0,2047,1)*dx)+dx/2-180. ; note -180. added by sjm to align
lon!0 = "lon"
lon&lon = lon
lon@units = "degrees_east"
;***************************************
; fill out the netCDF data model
;***************************************
sst!0 = "lat" ; name dimensions
sst!1 = "lon" ; ditto
sst = sst(::-1,:) ; reverse lat orientation
sst@long_name = "NAVO MCSST" ; assign long_name
sst@units = "deg C" ; assign units
sst&lat = lat ; assign lat cv
sst&lon = lon ; assign lon cv
sst@_FillValue = -999. ; assign missing value
;***************************************
; get year and day from filename
;***************************************
res = True ; plot mods desired
title = stringtochar(fname) ; parse file name to get date
year = title(0:3)
jday = title(4:6)
res@gsnCenterString = year+" "+jday ; create center string
;***************************************
; create plot
;***************************************
wks = gsn_open_wks("ps","mcsst") ; open workstation (plot destination)
gsn_define_colormap(wks,"BlGrYeOrReVi200") ; choose colormap
;
; This will not be necessary in V6.1.0 and later. Named colors can
; be used without having to first add them to the color map.
;
d = NhlNewColor(wks,0.8,0.8,0.8) ; add gray to colormap
res@cnFillOn = True ; turn on color
res@gsnSpreadColors = True ; use full range of colormap
res@gsnSpreadColorStart = 2 ; start at color 2
res@gsnSpreadColorEnd = -3 ; don't use added gray
res@cnLinesOn = False ; no contour lines
res@cnFillDrawOrder = "PreDraw" ; draw contours before continents
res@gsnMaximize = True ; maximize plot
; For a grid this size, it is better to use raster mode. It will be
; significantly faster, and will not go over NCL's 16mb default plot size.
res@cnFillMode = "RasterFill" ; turn on raster mode
plot = gsn_csm_contour_map_ce(wks,sst,res) ; contour the variable
end

3
samples/NCL/primero.ncl Normal file
View File

@@ -0,0 +1,3 @@
val=102
a=val/4.
print(a)

172
samples/NCL/topo_9.ncl Normal file
View File

@@ -0,0 +1,172 @@
;----------------------------------------------------------------------
; topo_9.ncl
;
; Concepts illustrated:
; - Recreating a jpeg topographic image as an NCL map object
; - Zooming in on a jpeg image
; - Drawing a box around an area of interest on a map
; - Attaching polylines to a map
; - Using "overlay" to overlay multiple contour plots
; - Using more than 256 colors per frame
; - Using functions for cleaner code
;----------------------------------------------------------------------
; NOTE: This example will only work with NCL V6.1.0 and later.
;
; This script recreates a JPEG image that was converted to a NetCDF
; file with color separated bands using the open source tool
; "gdal_translate":
;
; gdal_translate -ot Int16 -of netCDF EarthMap_2500x1250.jpg \
; EarthMap_2500x1250.nc
;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;----------------------------------------------------------------------
; This function imports a JPEG image that's on the whole globe,
; and recreates it as an NCL map object that is zoomed in on the
; southern tip of Africa.
;----------------------------------------------------------------------
undef("recreate_jpeg_image")
function recreate_jpeg_image(wks,minlat,maxlat,minlon,maxlon)
begin
orig_jpg_filename = "EarthMap_2500x1250.jpg"
nc_filename = "EarthMap_2500x1250.nc"
;--You could use a system call to do the NetCDF conversion
; cmd = "gdal_translate -ot Int16 -of netCDF " + jpeg_filename + \
; " " + nc_filename)
; system(cmd)
;---Read the three bands of data
f = addfile(nc_filename,"r")
Band1 = where(f->Band1.gt.255, 255, f->Band1) ; red channel
Band2 = where(f->Band2.gt.255, 255, f->Band2) ; green channel
Band3 = where(f->Band3.gt.255, 255, f->Band3) ; blue channel
band_dims = dimsizes(Band3)
nlat = band_dims(0)
nlon = band_dims(1)
print("dimensions of image = " + nlat + " x " + nlon)
;
; Add lat/lon data so we can overlay on a map, and/or
; overlay contours. We know the image is global,
; cylindrical equidistant, and centered about lon=0.
;
lat = fspan( -90, 90,nlat)
lon = fspan(-180,180,nlon)
lat@units = "degrees_north"
lon@units = "degrees_east"
Band1!0 = "lat"
Band1!1 = "lon"
Band2!0 = "lat"
Band2!1 = "lon"
Band3!0 = "lat"
Band3!1 = "lon"
Band1&lat = lat
Band1&lon = lon
Band2&lat = lat
Band2&lon = lon
Band3&lat = lat
Band3&lon = lon
res = True
res@gsnMaximize = True
res@gsnFrame = False ; Don't draw or advance
res@gsnDraw = False ; frame yet.
res@cnFillOn = True
res@cnFillMode = "RasterFill" ; Raster fill can be faster
res@cnLevelSelectionMode = "EqualSpacedLevels"
res@cnMaxLevelCount = 254
res@cnFillBackgroundColor = (/ 1., 1., 1., 1./)
res@cnLinesOn = False ; Turn off contour lines .
res@cnLineLabelsOn = False ; Turn off contour labels
res@cnInfoLabelOn = False ; Turn off info label
res@lbLabelBarOn = False ; Turn off labelbar
res@gsnRightString = "" ; Turn off subtitles
res@gsnLeftString = ""
res@pmTickMarkDisplayMode = "Always"
;---Construct RGBA colormaps...
ramp = fspan(0., 1., 255)
reds = new((/255, 4/), float)
greens = new((/255, 4/), float)
blues = new((/255, 4/), float)
reds = 0
greens = 0
blues = 0
reds(:,0) = ramp
greens(:,1) = ramp
blues(:,2) = ramp
; The red contour map is plotted fully opaque; the green and blue
; are plotted completely transparent. When overlain, the colors
; combine (rather magically).
reds(:,3) = 1.
greens(:,3) = 0
blues(:,3) = 0
res@cnFillColors = greens
greenMap = gsn_csm_contour(wks, Band2, res)
res@cnFillColors = blues
blueMap = gsn_csm_contour(wks, Band3, res)
;---This will be our base, so make it a map plot.
res@cnFillColors = reds
res@gsnAddCyclic = False
res@mpFillOn = False
;---Zoom in on area of interest
res@mpMinLatF = minlat
res@mpMaxLatF = maxlat
res@mpMinLonF = minlon
res@mpMaxLonF = maxlon
redMap = gsn_csm_contour_map(wks, Band1, res)
;---Overlay everything to create the topo map
overlay(redMap, greenMap)
overlay(redMap, blueMap)
return(redMap)
end
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
;---Recreating jpeg images only works for X11 and PNG.
wks = gsn_open_wks("png","topo")
;---Southern part of Africa
minlat = -40
maxlat = 5
minlon = 10
maxlon = 40
map = recreate_jpeg_image(wks,minlat,maxlat,minlon,maxlon)
;---Overlay a red box
lonbox = (/ 15, 35, 35, 15, 15/)
latbox = (/-30,-30,-10,-10,-30/)
lnres = True
lnres@gsLineColor = "red" ; red box
lnres@gsLineThicknessF = 4.0 ; make box thicker
box = gsn_add_polyline(wks,map,lonbox,latbox,lnres)
draw(map) ; Drawing the map will draw the red box
frame(wks)
end

120
samples/NCL/traj_3.ncl Normal file
View File

@@ -0,0 +1,120 @@
;*************************************************
; traj_3.ncl
;*************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
external TRAJ "./particle.so"
;*************************************************
begin
path = "./data.asc"
data = asciiread(path,(/500,6/),"float")
;*************************************************
; some parameters
;*************************************************
np = 1
nq = 500
ncor= 8
xrot = new((/np,nq/),float)
yrot = new((/np,nq/),float)
xaxis = new(ncor,float)
yaxis = new(ncor,float)
;**************************************************
; convert data into rotated format
;**************************************************
TRAJ::particle(path,xrot,yrot,nq,np,xaxis,yaxis,ncor)
;**************************************************
; create plot
;**************************************************
wks = gsn_open_wks("ps","traj") ; Open an ps file
xyres = True
xyres@gsnFrame = False ; don't advance the frame
xyres@gsnDraw = False ; don't draw indivdual plots
xyres@tmXTBorderOn = False ; don't draw top axis
xyres@tmXBBorderOn = False ; don't draw bottom axis
xyres@tmYRBorderOn = False ; don't draw right axis
xyres@tmYLBorderOn = False ; don't draw left axis
xyres@tmXTOn = False ; don't draw top-axis tick marks
xyres@tmXBOn = False ; don't draw bottom-axis tick marks
xyres@tmYROn = False ; don't draw right-axis tick marks
xyres@tmYLOn = False ; don't draw left-axis tick marks
xyres@xyLineColors = (/"red"/) ; set the line color to red
xyres@xyLineThicknessF = 4.0 ; 4 times the line thickness
xyres@trXMaxF = 15000 ; choose range of axis even though
xyres@trXMinF = -10000 ; we don't see them
xyres@trYMaxF = 1000
xyres@trYMinF = -1000
plot = gsn_xy(wks,xrot,yrot,xyres) ; Draw trajectory
;**********************************************
; create arrays needed for the bounding box
;**********************************************
a1 = new(5,float)
b1 = new(5,float)
a2 = new(5,float)
b2 = new(5,float)
a3 = new(2,float)
b3 = new(2,float)
a4 = new(2,float)
b4 = new(2,float)
a5 = new(2,float)
b5 = new(2,float)
a6 = new(2,float)
b6 = new(2,float)
a0 = new(2,float)
b0 = new(2,float)
;**********************************************
; determine values of each bounding line from information
; returned from particle.f
;**********************************************
a1(0:3) = xaxis(:3)
b1(0:3) = yaxis(:3)
a1(4) = xaxis(0)
b1(4) = yaxis(0)
a2(0:3) = xaxis(4:)
b2(0:3) = yaxis(4:)
a2(4) = xaxis(4)
b2(4) = yaxis(4)
a3 = xaxis(0:4:4)
b3 = yaxis(0:4:4)
a4 = xaxis(1:5:4)
b4 = yaxis(1:5:4)
a5 = xaxis(2:6:4)
b5 = yaxis(2:6:4)
a6 = xaxis(3:7:4)
b6 = yaxis(3:7:4)
a0(0) = xaxis(3)
b0(0) = yaxis(3)
a0(1) = xrot(0,0)
b0(1) = yrot(0,0)
;***************************************************************
; create bounding box by drawing multiple xy plots on top of
; each other. each with their individual axis turned off.
;***************************************************************
xyres@xyLineColors = (/"black"/) ; line color
xyres@xyLineThicknessF = 1.0 ; regular line thickness
bottom = gsn_xy(wks,a1,b1,xyres) ; Draw the bottom bounding box.
top = gsn_xy(wks,a2,b2,xyres) ; Draw the top bounding box.
side1 = gsn_xy(wks,a3,b3,xyres) ; Draw a side line.
side2 = gsn_xy(wks,a4,b4,xyres) ; Draw a side line.
side3 = gsn_xy(wks,a5,b5,xyres) ; Draw a side line.
side4 = gsn_xy(wks,a6,b6,xyres) ; Draw a side line.
;***************************************************************
; now draw a large brown line to represent the chimney
;***************************************************************
xyres@xyLineColors = (/"brown"/) ; chimney color
xyres@xyLineThicknessF = 9.0 ; thick line
xyres@tiMainString = "Pollutant Trajectory in a 3D Volume"
chimney = gsn_xy(wks,a0,b0,xyres) ; Draw the chimney.
draw(wks)
frame(wks)
end

167
samples/NCL/tsdiagram_1.ncl Normal file
View File

@@ -0,0 +1,167 @@
; Read potential temp (TEMP), salinity (SALT)
; Compute potential density (PD) for specified range PD(t,s)
; (use ncl function based on Yeager's algorithm for rho computation)
; Assumes annual and zonally avgeraged input data set (i.e, one time slice)
; Used K.Lindsay's "za" for zonal avg -- already binned into basins
; Plots temp vs salt (scatter plot), pd overlay
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
begin
; ================================> ; PARAMETERS
case = "PHC2_gx1v3"
ocnfile = "za_PHC2_T_S_gx1v3.nc"
depth_min = 14895.82 ; in cm, depth of first layer to be included
depth_max = 537499.9
;
; plot limits
;
smincn = 32.5
smaxcn = 37.0
tmincn = -2.
tmaxcn = 22.
;
; Choose basin index
;
; 0 = global 1 = southern ocean 2 = pacific 3 = indian 6 = atlantic
; 8 = labrador 9 = GIN 10 = arctic
;
bi = 2
;=====> basin check
if(bi.lt.0.or.bi.gt.10) then
print("basin index "+ bi + " not supported")
exit
end if
if(bi.eq.0) then
basin = "Global"
blab = "global"
end if
if(bi.eq.1) then
basin = "Southern Ocean"
blab = "so"
end if
if(bi.eq.2) then
basin = "Pacific Ocean"
blab = "pacific"
end if
if(bi.eq.3) then
basin = "Indian Ocean"
blab = "indian"
end if
if(bi.eq.6) then
basin = "Atlantic Ocean"
blab = "atlanticn"
end if
if(bi.eq.8) then
basin = "Labrador Sea"
blab = "lab"
end if
if(bi.eq.9) then
basin = "GIN Sea"
blab = "gin"
end if
if(bi.eq.10) then
basin = "Arctic Ocean"
blab = "arctic"
end if
;=====> initial resource settings
wks = gsn_open_wks("ps","tsdiagram") ; Open a Postscript file
;===== data
focn = addfile(ocnfile, "r")
salt = focn->SALT(0,:,{depth_min:depth_max},:) ;(basins, z_t, lat_t)
temp = focn->TEMP(0,:,{depth_min:depth_max},:)
;====section out choice basin
temp_ba = temp(bi,:,:)
salt_ba = salt(bi,:,:)
;===== put into scatter array format
tdata_ba = ndtooned(temp_ba)
sdata_ba = ndtooned(salt_ba)
ydata = tdata_ba
xdata = sdata_ba
;============== compute potenial density (PD), using rho_mwjf
;
; for potential density, depth = 0. (i.e. density as if brought to surface)
;
;===========================================================================
; WARNING: T-S diagrams use POTENTIAL DENSITY... if set depth to something
; other then 0, then you will be plotting density contours computed for the
; specified depth layer.
;===========================================================================
depth = 0. ;in meters
tspan = fspan(tmincn,tmaxcn,51)
sspan = fspan(smincn,smaxcn,51)
; the more points the better... using Yeager's numbers
t_range = conform_dims((/51,51/),tspan,0)
s_range = conform_dims((/51,51/),sspan,1)
pd = rho_mwjf(t_range,s_range,depth)
pd!0 = "temp"
pd!1 = "salt"
pd&temp = tspan
pd&salt = sspan
pd = 1000.*(pd-1.) ; Put into kg/m3 pot den units
; printVarSummary(pd)
; printVarInfo(pd,"rho_mwjf")
;=================Graphics
;--- scatter plot
res = True
res@gsnMaximize = True
res@gsnDraw = False
res@gsnFrame = False
res@xyMarkLineModes = "Markers"
res@xyMarkers = 16
res@xyMarkerColors = "black"
res@pmLegendDisplayMode = "Never"
res@txFontHeightF = 0.01
res@tiMainString = case + " ANN AVG: T-S Diagram"
res@tiXAxisString = salt@units
res@tiXAxisFontHeightF = 0.02
res@tiYAxisString = temp@units
res@tiYAxisFontHeightF = 0.02
res@trXMinF = smincn
res@trXMaxF = smaxcn
res@trYMinF = tmincn
res@trYMaxF = tmaxcn
res@gsnRightString = depth_min/100. + "-"+depth_max/100. +"m"
res@gsnLeftString = basin
plot = gsn_csm_xy(wks,xdata,ydata,res)
;----- pd overlay
resov = True
resov@gsnDraw = False
resov@gsnFrame = False
resov@cnLevelSelectionMode = "AutomaticLevels"
resov@cnInfoLabelOn = "False"
resov@cnLineLabelPlacementMode = "Constant"
resov@cnLineLabelFontHeightF = ".02"
plotpd = gsn_csm_contour(wks,pd,resov)
overlay(plot,plotpd)
draw(plot)
frame(wks)
end

141
samples/NCL/unique_9.ncl Normal file
View File

@@ -0,0 +1,141 @@
;************************************
; unique_9.ncl
;
; Concepts illustrated:
; - Drawing raster contours over a map
; - Creating a topography plot using raster contours
; - Reading data from binary files
; - Manually creating lat/lon coordinate arrays
; - Customizing a labelbar for a contour plot
;************************************
; This example generates a topo map over
; the area of Trinidad, Colorado.
;************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin
wks = gsn_open_wks("ps","unique")
;----------------- read the west binary data -------------------------
binfile = "trinidad-w.bin"
quad_name = fbinrecread(binfile,0,60,"character")
map_cornersW = fbinrecread(binfile,1,4,"double")
lonW = fbinrecread(binfile,2,(/1201/),"double")
latW = fbinrecread(binfile,3,(/1201/),"double")
minmax_elevW = fbinrecread(binfile,4,2,"double")
tmpW = fbinrecread(binfile,5,(/1201,1201/),"integer")
;----------------- read the east binary data -------------------------
binfile = "trinidad-e.bin"
quad_name = fbinrecread(binfile,0,60,"character")
map_cornersE = fbinrecread(binfile,1,4,"double")
lonE = fbinrecread(binfile,2,(/1201/),"double")
latE = fbinrecread(binfile,3,(/1201/),"double")
minmax_elevE = fbinrecread(binfile,4,2,"double")
tmpE = fbinrecread(binfile,5,(/1201,1201/),"integer")
;----------------------------------------------------------------------
min_elev = min((/minmax_elevW(0),minmax_elevE(0)/))*3.28
max_elev = max((/minmax_elevW(1),minmax_elevE(1)/))*3.28
lat = new(1201,"double")
lat = latW
lat!0 = "lat"
lat&lat = latW ; same as latE
lat@long_name = "latitude"
lat@units = "degrees_north"
lon = new(2401,"double")
lon(0:1200) = lonW
lon(1201:2400) = lonE(1:1200)
lon!0 = "lon"
lon&lon = lon
lon@long_name = "longitude"
lon@units = "degrees_east"
data = new((/1201,2401/),"float") ; (lat,lon)
data!0 = "lat"
data&lat = lat
data!1 = "lon"
data&lon = lon
data(:,0:1200) = (/tmpW*3.28/) ; convert to feet
data(:,1201:2400) = (/tmpE(:,1:1200)*3.28/) ; convert to feet
;-------------------------------------------------------------
;
; Define colormap.
;
cmap = (/(/1.00, 1.00, 1.00/),(/0.00, 0.00, 0.00/), \
(/0.51, 0.13, 0.94/),(/0.00, 0.00, 0.59/), \
(/0.00, 0.00, 0.80/),(/0.25, 0.41, 0.88/), \
(/0.12, 0.56, 1.00/),(/0.00, 0.75, 1.00/), \
(/0.63, 0.82, 1.00/),(/0.82, 0.96, 1.00/), \
(/1.00, 1.00, 0.78/),(/1.00, 0.88, 0.20/), \
(/1.00, 0.67, 0.00/),(/1.00, 0.43, 0.00/), \
(/1.00, 0.00, 0.00/),(/0.78, 0.00, 0.00/), \
(/0.63, 0.14, 0.14/),(/1.00, 0.41, 0.70/)/)
gsn_define_colormap(wks,cmap)
res = True
res@gsnMaximize = True
res@gsnAddCyclic = False
; map plot resources
res@mpFillOn = False
res@mpLimitMode = "Corners"
res@mpDataBaseVersion = "Ncarg4_1"
res@mpOutlineBoundarySets = "AllBoundaries"
res@mpLeftCornerLonF = map_cornersW(0)
res@mpLeftCornerLatF = map_cornersW(1)
res@mpRightCornerLonF = map_cornersE(2)
res@mpRightCornerLatF = map_cornersE(3)
; contour resources
res@cnFillOn = True
res@cnLinesOn = False
res@cnFillMode = "RasterFill"
res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels = (/ 5000., 6000., 7000., 8000., 8500., 9000., \
9500.,10000.,10500.,11000.,11500.,12000., \
12500.,13000.,13500./)
; tickmark resources
res@pmTickMarkDisplayMode = "Always"
res@tmXBLabelFontHeightF = 0.010
; labelbar resources
res@pmLabelBarWidthF = 0.60
res@txFontHeightF = 0.012
res@lbTitleString = "elevation above mean sea level (feet)"
res@lbTitleFontHeightF = 0.012
res@lbLabelFontHeightF = 0.008
res@lbTitleOffsetF = -0.27
res@lbBoxMinorExtentF = 0.15
res@pmLabelBarOrthogonalPosF = -.05
; title resources
res@tiMainString = "USGS DEM TRINIDAD (1 x 2 degrees)"
res@tiMainOffsetYF = -0.02 ; Move title down towards graphic.
res@tiMainFontHeightF = 0.015
res@gsnLeftString = "Min Elevation: "+min_elev
res@gsnRightString = "Max Elevation: "+max_elev
res@gsnCenterString = "Scale 1:250,000"
plot = gsn_csm_contour_map(wks,data,res)
end

131
samples/NCL/viewport_4.ncl Normal file
View File

@@ -0,0 +1,131 @@
; ***********************************************
; viewport_4.ncl
;
; Concepts illustrated:
; - Drawing an XY plot with multiple curves
; - Using drawNDCGrid to draw a nicely labeled NDC grid
; - Changing the size/shape of an XY plot using viewport resources
; - Drawing two XY plots on the same page using viewport resources
; - Drawing polylines, polymarkers, and text in NDC space
; - Using "getvalues" to retrieve resource values
; - Maximizing plots after they've been created
; ***********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;********************************************************************
; Draw a box around the viewport of the given object..
;********************************************************************
procedure draw_vp_box(wks,plot)
local vpx, vpy, vpw, vph, xbox, ybox, lnres, mkres, txres
begin
; Retrieve the viewport values of the drawable object.
getvalues plot
"vpXF" : vpx
"vpYF" : vpy
"vpWidthF" : vpw
"vpHeightF" : vph
end getvalues
; Set up some marker resources.
mkres = True
mkres@gsMarkerIndex = 16 ; filled dot
mkres@gsMarkerSizeF = 0.02 ; larger than default
mkres@gsMarkerColor = "Red"
; Draw a single marker at the vpXF/vpYF location.
gsn_polymarker_ndc(wks,vpx,vpy,mkres)
; Set up some text resources.
txres = True
txres@txJust = "BottomLeft"
txres@txFontHeightF = 0.018
txres@txFontColor = "Blue"
txres@txBackgroundFillColor = "white"
gsn_text_ndc(wks,"(vpXF="+vpx+", vpYF="+vpy+")",vpx,vpy+0.02,txres)
; Set up some line resources.
lnres = True
lnres@gsLineColor = "Red" ; line color
lnres@gsLineThicknessF = 2.0 ; 3.5 times as thick
; Draw lines indicating the width and height
xline = (/vpx, vpx+vpw/)
yline = (/vpy-0.05,vpy-0.05/)
gsn_polyline_ndc(wks,xline,yline,lnres)
xline = (/vpx+0.05,vpx+0.05/)
yline = (/vpy,vpy-vph/)
gsn_polyline_ndc(wks,xline,yline,lnres)
txres@txJust = "CenterCenter"
gsn_text_ndc(wks,"vpWidthF = " + vpw,vpx+vpw/2.,vpy-0.05,txres)
txres@txAngleF = 90.
gsn_text_ndc(wks,"vpHeightF = " + vph,vpx+0.05,vpy-vph/2.,txres)
end
;********************************************************************
; Main code
;********************************************************************
begin
;************************************************
; read in data
;************************************************
f = addfile ("$NCARG_ROOT/lib/ncarg/data/cdf/uv300.nc","r")
u = f->U ; get u data
;************************************************
; plotting parameters
;************************************************
wks = gsn_open_wks ("ps","viewport") ; open workstation
res = True ; plot mods desired
res@gsnFrame = False ; don't advance frame yet
res@vpWidthF = 0.8 ; set width and height
res@vpHeightF = 0.3
; First plot
res@tiMainString = "Plot 1"
res@vpXF = 0.15
res@vpYF = 0.9 ; Higher on the page
plot1 = gsn_csm_xy (wks,u&lat,u(0,:,{82}),res) ; create plot
; Second plot
res@tiMainString = "Plot 2"
res@vpXF = 0.15 ; Same X location as first plot
res@vpYF = 0.4 ; Lower on the page
plot2 = gsn_csm_xy (wks,u&lat,u(0,:,{3}),res) ; create plot
; Advance the frame
frame(wks)
; Now draw the two plots with illustrations.
drawNDCGrid(wks) ; Draw helpful grid lines showing NDC square.
draw(plot1) ; Draw the two plots
draw(plot2)
draw_vp_box(wks,plot1) ; Draw boxes around the two viewports.
draw_vp_box(wks,plot2)
frame(wks) ; Advance the frame.
;
; Uncomment the next two lines if you want to maximize these plots for
; PS or PDF output.
;
; psres = True
; maximize_output(wks,psres) ; calls draw and frame for you
end

120
samples/NCL/weather_sym_6.ncl Normal file
View File

@@ -0,0 +1,120 @@
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; Example of plotting station model data over a map
; illustrating how the wind barb directions are adjusted
; for the map projection.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; City names.
;
cities = (/ "NCAR", "Seattle", "San Francisco", \
"Los Angeles", "Billings", "El Paso", \
"Houston", "Kansas City", "Minneapolis", \
"Chicago", "Detroit", "Atlanta", \
"Miami", "New York", "Eugene", \
"Boise", "Salt Lake", "Phoenix", \
"Albuquerque", "Bismarck", "Tulsa", \
"Dallas", "Little Rock", "Lexington", \
"Charlotte", "Norfolk", "Bangor" \
/)
city_lats = (/ 40.0, 47.6, 37.8, \
34.1, 45.8, 31.8, \
29.8, 39.1, 45.0, \
41.9, 42.3, 33.8, \
25.8, 40.8, 44.1, \
43.6, 40.7, 33.5, \
35.1, 46.7, 36.0, \
32.8, 34.7, 38.1, \
35.2, 36.8, 44.8 \
/)
city_lons = (/ -105.0, -122.3, -122.4, \
-118.3, -108.5, -106.5, \
-095.3, -094.1, -093.8, \
-087.6, -083.1, -084.4, \
-080.2, -074.0, -123.1, \
-116.2, -111.9, -112.1, \
-106.6, -100.8, -096.0, \
-096.8, -092.3, -084.1, \
-080.8, -076.3, -068.8 \
/)
;
; Station model data for the 27 cities.
;
imdat = (/"11000000751126021360300004955054054600007757087712", \
"11103100011104021080300004959055050600517043080369", \
"11206200031102021040300004963056046601517084081470", \
"11309300061000021020300004967057042602017125082581", \
"11412400091002021010300004971058038602517166083592", \
"11515500121004020000300004975050034603017207084703", \
"11618600151006020030300004979051030603507248085814", \
"11721700181008020050300004983052026604007289086925", \
"11824800211009020070300004987053022604507323087036", \
"11927900241011020110300004991054018605017364088147", \
"11030000271013020130300004995055014605517405089258", \
"11133100301015020170300004999056010606017446080369", \
"11236200331017020200300004000057006606517487081470", \
"11339300361019020230300004004058002607017528082581", \
"11442400391021020250300004008050000607517569083692", \
"11545500421023020270300004012051040608017603084703", \
"11648600451025020290300004017052008608517644085814", \
"11751700481027020310300004021053012609017685086925", \
"11854800511029020330300004025054016609507726087036", \
"11958900541031020360300004029055018610007767088147", \
"11060000571033020380300004033056030610507808089258", \
"11163100601035020410300004037057034611007849080369", \
"11266200631037020430300004041058043611507883081470", \
"11369300661039020470300004045050041612007924082581", \
"11472400691041020500300004048051025612507965083692", \
"11575500721043020530300004051052022613507996084703", \
"11678600751048021580300004055053013614007337085814" \
/)
;
; Define a color map and open a workstation.
;
cmap = (/ \
(/ 1., 1., 1. /), \ ; color index 0 - white
(/ 0., 0., 0. /) \ ; color index 1 - black
/)
wks = gsn_open_wks("ps","weather_sym")
gsn_define_colormap(wks,cmap)
;
; Draw a world map.
;
mpres = True
mpres@gsnFrame = False
mpres@mpSatelliteDistF = 1.3
mpres@mpOutlineBoundarySets = "USStates"
mpres@mpCenterLatF = 40.
mpres@mpCenterLonF = -97.
mpres@mpCenterRotF = 35.
map = gsn_map(wks,"Satellite",mpres)
;
; Scale the station model plot (all aspects of the station
; model plots are scaled as per the size of the wind barb).
;
wmsetp("wbs",0.018)
;
; In the middle of Nebraska, draw a wind barb for a north wind
; with a magnitude of 15 knots.
;
wmbarbmap(wks,42.,-99.,0.,-15.)
;
; Draw the station model data at the selected cities. The call
; to wmsetp informs wmstnm that the wind barbs will be drawn over
; a map. To illustrate the adjustment for plotting the model
; data over a map, all winds are from the north.
;
wmsetp("ezf",1)
wmstnm(wks,city_lats,city_lons,imdat)
frame(wks)
end

151
samples/NCL/xy_29.ncl Normal file
View File

@@ -0,0 +1,151 @@
; xy_29.ncl
;
; Concepts illustrated:
; - Reading data from an ASCII file with headers
; - Creating a separate procedure to create a specific plot
; - Attaching polymarkers to an XY plot
;
; This script was originally from Dr. Birgit Hassler (NOAA)
;****************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************
; Plot Procedure
;************************************************
procedure plotTCOPolym(pltName[1]:string, pltType[1]:string, filName[1]:string \
,xTitle[1]:string , yTitle[1]:string \
,year[*]:numeric, y[*]:numeric)
local wks, res, ntim, gsres, MarkerCol, OldYear, i, xmarker, ymarker
begin
wks = gsn_open_wks(pltType,pltName)
gsn_define_colormap(wks,"default")
res = True
res@gsnMaximize = True ; make "ps", "eps", "pdf" large
res@vpHeightF = 0.5 ; change aspect ratio of plot
res@vpWidthF = 0.75
res@vpXF = 0.15 ; start plot at x ndc coord
res@tiXAxisString = xTitle
res@tiYAxisString = yTitle
res@tiMainString = filName
ntim = dimsizes(year)
res@trXMinF = year(0)-1
res@trXMaxF = year(ntim-1)+1
res@gsnDraw = False
res@gsnFrame = False
res@xyMarkLineMode = "markers"
res@xyMarker = 16
res@xyMarkerColor = "Background"
plot = gsn_csm_xy (wks,year,y,res) ; create plot frame ork
; add different color polymarkers for each year
gsres = True
MarkerCol = 2
OldYear = year(0)
do i=0,ntim-1
xmarker = year(i)
ymarker = y(i)
if (i.gt.0) then
if (year(i).gt.OldYear) then
MarkerCol = MarkerCol+1
end if
OldYear = year(i)
end if
gsres@gsMarkerColor = MarkerCol
gsres@gsMarkerIndex = 16
;gsres@gsMarkerSizeF = 15.0
; add (attach) polymarkers to existing plot object
plot@$unique_string("dum")$ = gsn_add_polymarker(wks,plot,xmarker,ymarker,gsres)
end do
draw(plot)
frame(wks)
end
;***********************************************************
; MAIN
;***********************************************************
pltType = "ps" ; "ps", "eps", "png", "x11"
; read multiple ascii file names
;;fili = "Southpole_TCOTimeSeries_11.dat"
diri = "./"
fili = systemfunc("cd "+diri+" ; ls *TCOT*dat")
print(fili)
nfil = dimsizes(fili)
nhead= 4 ; number of header lines on ascii file(s)
ncol = 4 ; year, month, day, O3
do nf=0,nfil-1
sfx = get_file_suffix(fili(nf), 0) ; sfx = ".dat"
filx = sfx@fBase ; filx= "Southpole_TCOTimeSeries_11"
; read ascii files
data = readAsciiTable(diri+fili(nf), ncol, "float", nhead)
dimd = dimsizes(data)
ntim = dimd(0) ; # rows
year = toint( data(:,0) ) ; user decision ... convert to integer
mon = toint( data(:,1) )
day = toint( data(:,2) )
hour = new (ntim, "integer", "No_FillValue")
mn = new (ntim, "integer", "No_FillValue")
sec = new (ntim, "double" , "No_FillValue")
hour = 0
mn = 0
sec = 0d0
; create COARDS/udunits time variable
;;tunits = "days since 1900-01-01 00:00:0.0"
tunits = "days since "+year(0)+"-"+mon(0)+"-"+day(0)+" 00:00:0.0"
time = cd_inv_calendar(year,mon,day,hour,mn,sec,tunits, 0)
time!0 = "time"
time&time = time
;printVarSummary(time)
; create a Gregorin 'date' variable
date = year*10000 + mon*100 + day
date!0 = "time"
date@units = "yyyymmdd"
date&time = time
;printVarSummary(date)
O3 = data(:,3)
O3@long_name = "total column ozone"
O3@units = "DU"
O3!0 = "time"
O3&time = time
;printVarSummary(O3)
;print(" ")
;print(date+" "+time+" "+O3)
; plot
yTitle = O3@long_name
year@long_name = "YEAR"
plotTCOPolym (filx, pltType, fili(nf), year@long_name, yTitle, year, O3)
delete(time) ; delete ... size (# rows) may change in the next file
delete(date)
delete(year)
delete(mon )
delete(day )
delete(mn )
delete(sec )
delete(O3 )
delete(data)
end do