Re: Adding a shapefile to multiple plots

From: Mary Haley <haley_at_nyahnyahspammersnyahnyah>
Date: Tue Apr 24 2012 - 09:54:31 MDT

Hi Nkese,

I'm confused by something, because you reference both "plot" and "plots" in your code below.

I'm going to go with "plots" and guess that the problem is that you are passing in an array of plots to attach_shapefile_polyline, when it is expecting only one plot.

First: since you are attaching the polylines to "plot" as an attribute before you leave attach_shapefile_polyline,
you should change this to a procedure instead:

> undef("attach_shapefile_polyline")
> procedure attach_shapefile_polyline(f:file,wks,plot,line_color)

Then,, instead of this:

> plot = attach_shapefile_polyline(f,wks,plot,"black")
> plot = attach_shapefile_polyline(f1,wks,plot,"black")
>
> draw(plots)
> frame(wks)

Try this:

attach_shapefile_polyline(f,wks,plot(nmo),"black")
attach_shapefile_polyline(f1,wks,plot(nmo),"black")

Also, I'm not sure why you have "draw(plots)" when you later have
a call to gsn_panel. I would remove the "draw(plots)" line, because is going
to draw all of your plots on one frame, one on top of the other.
You will also need to remove the "frame(wks)" call.

--Mary

On Apr 23, 2012, at 4:10 PM, Nkese Mc Shine wrote:

> 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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Apr 24 09:54:37 2012

This archive was generated by hypermail 2.1.8 : Mon Apr 30 2012 - 09:21:12 MDT