Dear All,
Below is my attached script. It reads in 25 years of monthly data and produces for each year 12 maps indicating each month. I would like to know how can I add shapefiles to all these plots. When I try I am getting an error such as :
fatal:Number of elements of dimension (0) of argument (0) is (12) in function (NhlAddPrimitive), expected (1) elements
fatal:Execute: Error occurred at or near line 4147 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl
fatal:Execute: Error occurred at or near line 4412 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl
fatal:Execute: Error occurred at or near line 59 in file practice.ncl
fatal:Execute: Error occurred at or near line 392 in file practice.ncl
;----------------------------------------------------
; Fixing Quality Flag
;------------------------------------------------
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"
;----------------------------------------------------------------------
; This procedure attaches polyline data from a shapefile
;----------------------------------------------------------------------
undef("attach_shapefile_polyline")
function attach_shapefile_polyline(f:file,wks,plot,line_color)
local lnres
begin
; gsn_define_colormap(wks,(/"white","black","tan","LightBlue",\
; "brown","yellow","navyblue","green"/))
lnres = True ; resources for polylines
lnres@gsLineThicknessF = 2.0 ; 2x as thick
lnres@gsLineColor = line_color
;
; Loop through files that we want to read geographic information from.
;
; If this loop is extremely slow, consider using gsn_polyline instead
; of gsn_add_polyline. This can have a significant effect. Remember
; that gsn_polyline is a procedure, not a function, and it draws the
; lines right when you call it, so you need to make sure your plot is
; already drawn to the frame.
;
;---Read data off the shapefile
segments = f->segments
geometry = f->geometry
segsDims = dimsizes(segments)
geomDims = dimsizes(geometry)
;---Read global attributes
geom_segIndex = f@geom_segIndex
geom_numSegs = f@geom_numSegs
segs_xyzIndex = f@segs_xyzIndex
segs_numPnts = f@segs_numPnts
numFeatures = geomDims(0)
;---Section to draw polylines on plot.
lon = f->x
lat = f->y
do i=0, numFeatures-1
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
;
; This call adds the line segment.
;
; Can use gsn_polyline here to make it faster.
;
dumstr = unique_string("primitive")
plot@$dumstr$ = gsn_add_polyline(wks, plot, lon(startPT:endPT), \
lat(startPT:endPT), lnres)
end do
end do
;---Clean up before we read in same variables again.
delete(lat)
delete(lon)
delete(segments)
delete(geometry)
;---We have to return plot so that attached lines are not lost.
return(plot)
end
;-----------------------------------------------------
;MAIN SCRIPT
;------------------------------------------------------
begin
; specify region boundaries
latS = 9.5
latN = 11.5
lonL = -62
lonR = -60
; directories and files
DIR = "/media/data1/nfs_mount1/Data/"
DN = "day" ; "day" or "night"
;dirs = "day1986" ; one year
dirs = DN+"[1-2][0-9]*" ; all years
print(dirs)
dir_year = systemfunc("cd "+DIR+" ; ls -d "+dirs)
print(dir_year)
diri = "./"+dirs+"/"
fils = systemfunc("cd "+DIR+" ; ls "+diri+"*.nc")
;print(fils)
fad = addfiles (DIR+fils,"r")
ListSetType(fad,"join")
; Read global lat/lon from 1st file; find region index
LAT = (/ fad[0]->lat /)
iLAT = ind(LAT.le.latN .and. LAT.ge.latS)
nStrt = iLAT(0) ; start index
nLast = iLAT(dimsizes(iLAT)-1) ; last index
LON = (/ fad[0]->lon /)
iLON = ind(LON.le.lonR .and. LON.ge.lonL)
mStrt = iLON(0) ; start index
mLast = iLON(dimsizes(iLON)-1) ; last index
;print("nStrt="+nStrt+" nLast="+nLast)
;print("mStrt="+mStrt+" mLast="+mLast)
; Read from just the desired region
sst = (/ short2flt( fad[:]->sst(:,nStrt:nLast,mStrt:mLast) ) /) ; return no meta data
; manually create 'good' meta data
sst_copy = sst
lat = LAT(nStrt:nLast)
lat!0 = "lat"
lat@units = "degrees_north"
lon = LON(mStrt:mLast)
lon!0 = "lon"
lon@units = "degrees_east"
nyear = dimsizes(dir_year)
if (DN.eq."day") then
yrStrt = toint( str_get_cols(dir_year(0),3,6) )
yrLast = toint( str_get_cols(dir_year(nyear-1),3,6) )
else
yrStrt = toint( str_get_cols(dir_year(0),5,8) )
yrLast = toint( str_get_cols(dir_year(nyear-1),5,8) )
end if
time = yyyymm_time(yrStrt, yrLast, "integer")
sst!0 = "time"
sst!1 = "lat"
sst!2 = "lon"
sst_copy!0 = "time"
sst_copy!1 = "lat"
sst_copy!2 = "lon"
sst&time= time
sst&lat = lat
sst&lon = lon
sst_copy&time = time
sst_copy&lat = lat
sst_copy&lon = lon
sst@long_name = "SST"
sst@units = "degC"
sst@_FillValue = 1e20
sst_copy@_FillValue = sst@_FillValue
;sst = where(sst.le.-3, sst@_FillValue, sst) ;fix this depending on the results obtained from this script
;sst = where(sst.le.-3 .and. sst.eq.0, sst@_FillValue, sst) ;fix this depending on the results obtained from this script
printVarSummary(sst)
printMinMax(sst, True)
;print(sst(:,:,{-61}))
;---------------------------------
;Standard deviation for sst
;---------------------------------
;stddev_sst = stddev(sst)
;print(stddev_sst)
sststdTime = dim_stddev_n(sst,0)
;print(sststdTime)
;---------------------------------
;Standard deviation for sst_copy
;---------------------------------
;stddev_sst_copy = stddev(sst_copy)
;print(stddev_sst_copy)
sst_copy_stdTime = dim_stddev_n(sst_copy,0)
;print(sst_copy_stdTime)
; directories and files for qual
DIRqual = "/media/data1/nfs_mount1/Qualityflag/"
DNqual = "qday"
;dirsqual="qday1985"
dirsqual= DNqual+"[1-2][0-9]*"
print(dirsqual)
dir_yearqual = systemfunc("cd "+DIRqual+" ; ls -d "+dirsqual)
print(dir_yearqual)
diriqual = "./"+dirsqual+"/"
filsqual = systemfunc("cd "+DIRqual+" ; ls "+diriqual+"*.nc")
;print(filsqual)
fadqual = addfiles(DIRqual+filsqual,"r")
ListSetType(fadqual,"join")
; Read global lat/lon from 1st file; find region index for qsst
LATQ = (/ fadqual[0]->lat /)
iLAT = ind(LAT.le.latN .and. LAT.ge.latS)
nStrt = iLAT(0) ; start index
nLast = iLAT(dimsizes(iLAT)-1) ; last index
LONQ = (/ fadqual[0]->lon /)
iLON = ind(LON.le.lonR .and. LON.ge.lonL)
mStrt = iLON(0) ; start index
mLast = iLON(dimsizes(iLON)-1) ; last index
;print("nStrt="+nStrt+" nLast="+nLast)
;print("mStrt="+mStrt+" mLast="+mLast)
; Read from just the desired region for qsst
qsst = (/ ( fadqual[:]->qual(:,nStrt:nLast,mStrt:mLast) ) /) ;return no meta data
;qsst = fadqual[:]->qual(nStrt:nLast,mStrt:mLast)
;qsst = fadqual[:]->qual({latS:latN},{lonL:lonR})
printVarSummary(qsst)
; dsp_Flag = qsst
latq = LATQ(nStrt:nLast)
latq!0 = "lat"
latq@units = "degrees_north"
lonq = LONQ(mStrt:mLast)
lonq!0 = "lon"
lonq@units = "degrees_east"
qnyear = dimsizes(dir_yearqual)
; print(nStrt)
; print(nLast)
; print(qnyear)
; aaa=str_get_cols(dir_yearqual(0),4,7)
; print(aaa)
; print(dir_yearqual(0))
if (DNqual.eq."qday") then
yrStrt = toint( str_get_cols(dir_yearqual(0),4,7) )
yrLast = toint( str_get_cols(dir_yearqual(qnyear-1),4,7) )
else
yrStrt = toint( str_get_cols(dir_yearqual(0),6,9) )
yrLast = toint( str_get_cols(dir_yearqual(qnyear-1),6,9) )
end if
qtime = yyyymm_time(yrStrt, yrLast, "integer")
;print(qtime)
qsst!0 = "qtime"
qsst!1 = "lat"
qsst!2 = "lon"
qsst&qtime= qtime
qsst&lat = lat
qsst&lon = lon
printVarSummary(qsst)
dimqsst = dimsizes(qsst)
nlat = dimqsst(1)
nlon = dimqsst(2)
;print(dimqsst)
; Checking through each cell, find the times at which quality flag is not good. Range 0-4
do i=1,nlat
do j=1,nlon
sst_copy = where(qsst.ge.0 .and. qsst.le.4, sst_copy@_FillValue, sst_copy)
;notgoodtimes = where(qsst(:,i,j).ge.0 .or. qsst(:,i,j).le.4)
;qsst(notgoodtimes,i,j) = sst@_FillValue
end do
end do
printVarSummary(sst_copy)
printMinMax(sst_copy, True)
;print(sst_copy(:,:,{-61}))
;creating a plot
;wks = gsn_open_wks("x11" ,"sst") ; open a ps file
wks = gsn_open_wks("pdf" ,"sst")
setvalues NhlGetWorkspaceObjectId()
"wsMaximumSize" : 300000000
end setvalues
;---Force colors 2-7 to be tan
tan_rgb = namedcolor2rgb("tan")
NhlSetColor(wks,2,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
NhlSetColor(wks,3,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
NhlSetColor(wks,4,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
NhlSetColor(wks,5,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
NhlSetColor(wks,6,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
NhlSetColor(wks,7,tan_rgb(0,0),tan_rgb(0,1),tan_rgb(0,2))
;---Draw the color map
; gsn_draw_colormap(wks)
gsn_define_colormap(wks,"rainbow") ; choose colormap
res = True ; plot mods desired
res@gsnMaximize = True ; make ps, eps, pdf large
res@gsnDraw = False ; don't draw plot
res@gsnFrame = False ; don't advance frame
res@cnFillOn = True ; turn on color fill
res@cnFillMode = "RasterFill" ; Raster Mode
res@cnLinesOn = False ; turn of contour lines
res@cnLineLabelsOn = False ; turn of contour line labels
res@cnLevelSelectionMode = "EqualSpacedLevels"
res@cnMinLevelValF = 26.5
res@cnMaxLevelValF = 35
res@cnLevelSpacingF = 0.5 ; contour spacing
res@gsnSpreadColors = True ; use full range of color map
res@lbOrientation = "Vertical";Vertical labelbar
;res@lbLabelStride = 4
res@lbLabelAutoStride = False ; let NCL choose spacing
res@lbLabelBarOn = False
;plot = gsn_csm_contour_map_ce(wks,sst, res) ; global
res@mpOutlineOn = False ; Turn off map outlines
res@mpOceanFillColor = 11
res@mpLandFillColor = 192
res@mpInlandWaterFillColor = 74
res@mpMinLatF = latS
res@mpMaxLatF = latN
res@mpMinLonF = lonL
res@mpMaxLonF = lonR
res@gsnAddCyclic = False ; not global
;resP@txString = " DAY SST"
; plot = gsn_csm_contour_map_ce(wks,sst_copy(0,{latS:latN},{lonL:lonR}),res) ;regional
;---Open two shapefiles.
; f = addfile("TTO_adm0.shp", "r")
; f1= addfile("VEN_adm0.shp", "r")
;---Attach polylines to map.
; plot = attach_shapefile_polyline(f,wks,plot,"black")
; plot = attach_shapefile_polyline(f1,wks,plot,"black")
; draw(plot)
; frame(wks)
plots = new(12,graphic)
plot = plots
months = (/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug",\
"Sep","Oct","Nov","Dec"/)
years = "" + ispan(1985,2009,1)
pres = True ; Set panel resources
pres@gsnMaximize = True ; Maximize plots in panel
pres@gsnPanelLabelBar = True ; Turn on panel labelbar
pres@lbLabelAutoStride = True ; Nice spacing for labels
pres@lbBoxLinesOn = False
do nyr=0,nyear-1
do nmo=0,11
imo = nyr*12 + nmo
res@gsnCenterString=months(nmo)
plots(nmo) = gsn_csm_contour_map_ce(wks,sst_copy(imo,{latS:latN},{lonL:lonR}),res)
end do
pres@txString = years(nyr)
;plots = gsn_csm_contour_map_ce(wks,sst_copy(imo,{latS:latN},{lonL:lonR}),res) ;regional
;---Open two shapefiles.
f = addfile("TTO_adm0.shp", "r")
f1= addfile("VEN_adm0.shp", "r")
;---Attach polylines to map.
plot = attach_shapefile_polyline(f,wks,plot,"black")
plot = attach_shapefile_polyline(f1,wks,plot,"black")
draw(plots)
frame(wks)
;create a panel of plots with 3 rows and 4 coloums
gsn_panel(wks,plots,(/3,4/),pres)
end do
end
Best regards,
Nkese
CONFIDENTIALITY: This email (including any attachments) may contain confidential, proprietary and/or privileged information. Any duplication, copying, distribution, dissemination, transmission, disclosure or use in any manner of this email (including any attachments) without the authorisation of the sender is strictly prohibited. If you receive this email (including any attachments) in error, please notify the sender and delete this email (including any attachments) from your system. Thank you.
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 23 16:10:30 2012
This archive was generated by hypermail 2.1.8 : Mon Apr 30 2012 - 09:21:12 MDT