Re: Masking & shapefiles

Hi Nkese,

It would help if you can include an image. I'm not sure what you mean by the shapefile scheme controlling the color scheme for the maps. You are drawing the shapefile outlines in black, so maybe what you think is the NCL map outline (which is drawn in black by default), is actually the shapefile outline?

As for your second question, you created a "sst_mask" variable, but then you never use it. Instead, you are plotting "sst" in your gsn_csm_contour_map call. Try changing this to plot sst_mask.

--Mary

On Feb 22, 2012, at 9:45 AM, Nkese Mc Shine wrote:

>
> 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

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Feb 22 13:34:26 2012