Dear All,
Attached is my script. It includes two shape-files, a masking script and latitude and longitude for a specific region. I have two problems:
1) The shape-files colour scheme controls the colour scheme for the maps instead of the graphical colour scheme for the maps.
2) I am trying to mask the ocean only but i am not sure where I'm putting the masking script in my script is in the correct position because I am not getting any graphical out puts.
Regards,
Nkese.
;---------------------------------------------------------
; REGION
;---------------------------------------------------------
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"
;----------------------------------------------------------------------
; 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,"rainbow")
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 = "/root/Desktop/NCLdatasets/Datasst1/"
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
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(0),5,8) )
end if
time = yyyymm_time(yrStrt, yrLast, "integer")
sst!0 = "time"
sst!1 = "lat"
sst!2 = "lon"
sst&time= time
sst&lat = lat
sst&lon = lon
sst@long_name = "SST"
sst@units = "degC"
sst@_FillValue = 1e20
sst = where(sst.le.-3, sst@_FillValue, sst) ;fix this depending on the results obtained from this script
;------------------------------------------------------
; End Required User Specifications
; Note: User may edit the mask below
;------------------------------------------------------
;---try to add mask to the data with a shapefile
nmrb = dimsizes(lon)
min_lat = min(lat)
max_lat = max(lat)
min_lon = min(lon)
max_lon = max(lon)
;
;---Get the approximate index values that contain the area of interest
;
;---This will make our gc_inout loop below go faster,beacuse we're not
; checking every single lat/lon point to see if it's within the area
; of interest.
;
ilt_mn = ind(min_lat.gt.lat)
ilt_mx = ind(max_lat.lt.lat)
iln_mn = ind(min_lon.gt.lon)
iln_mx = ind(max_lon.lt.lon)
ilt1 = ilt_mn(dimsizes(ilt_mn)-1) ;Start of lat box
iln1 = iln_mn(dimsizes(iln_mn)-1) ;Start of lon box
ilt2 = ilt_mx(0) ;End of lat box
iln2 = iln_mx(0) ;End of lon box
sst_mask = new(dimsizes(sst),typeof(sst),sst@_FillValue)
copy_VarCoords(sst,sst_mask)
;---Put data in the areas that we don't want masked.
do ilt=ilt1,ilt2
do iln=iln1,iln2
if(gc_inout(lat(ilt),lon(iln),lat,lon)) then
sst_mask(ilt,iln) = data(ilt,iln)
end if
end do
end do
printVarSummary(sst)
printMinMax(sst, True)
;print(sst(:,:,{-61}))
;;setfileoption("nc", "Format", "NetCDF4Classic") ; write variables > 2GB
;diro = "./"
;filo = DN+"."+yrStrt+"-"+yrLast+".nc"
;patho = diro+filo
;system("/bin/rm -f "+patho)
;ncdf = addfile(patho, "c")
;ncdf->SST = sst
;creating a plot
wks = gsn_open_wks("x11" ,"sst") ; open a ps file
setvalues NhlGetWorkspaceObjectId()
"wsMaximumSize" : 300000000
end setvalues
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 = "ManualLevels"
res@cnMinLevelValF = 15
res@cnMaxLevelValF = 38
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 = True ; let NCL choose spacing
res@lbLabelBarOn = True
;plot = gsn_csm_contour_map_ce(wks,sst, res) ; global
res@mpOutlineOn = False ; Turn off map outlines
res@mpOceanFillColor = 5
res@mpLandFillColor = 164
res@mpInlandWaterFillColor = 54
res@mpDataBaseVersion = "MediumRes" ; slightly better resolution
res@mpMinLatF = latS
res@mpMaxLatF = latN
res@mpMinLonF = lonL
res@mpMaxLonF = lonR
res@gsnAddCyclic = False ; not global
plot = gsn_csm_contour_map_ce(wks,sst(5,{latS:latN},{lonL:lonR}),res)
;plot(0) = gsn_csm_contour_map_ce(wks,sst(1,{latS:latN},{lonL:lonR}),res)
f1 = addfile("TTO_adm1.shp", "r")
f2 = addfile("VEN_adm0.shp", "r")
;---Attach polylines to map.
plot = attach_shapefile_polyline(f1,wks,plot,"black")
plot = attach_shapefile_polyline(f2,wks,plot,"black")
draw(plot)
frame(wks)
end
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 Wed Feb 22 09:45:45 2012
This archive was generated by hypermail 2.1.8 : Thu Feb 23 2012 - 10:01:53 MST