[ncl-talk] Mask using shapefile
Mary Haley
haley at ucar.edu
Wed Mar 8 16:36:07 MST 2017
Hi,
Can you provide the b.e11 NetCDF file? You can email me offline about this
if you want.
I do recognize the code you included because it's code that I wrote, but
it's hard to debug without the data file.
Also, did you look at the "shapefile_mask_data" function on the shapefiles
example page?
In particular, you could look at example shapefile_11.ncl:
http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex11
You can first try it with minimal options to see what happens:
opt = True
opt at debug = True
ts_mask = shapefile_mask_data(ts,India_shp,opt)
printVarSummary(ts_mask)
printMinMax(ts_mask,0)
I can't remember, but you might need to call this too:
copy_VarMeta(ts,ts_mask) ; Copy all metadata
--Mary
On Wed, Mar 8, 2017 at 12:59 PM, Ipshita Majhi <ipmajhi at alaska.edu> wrote:
> Dear NCL Users,
>
> I am trying to create a mask for India, using shapefile. I am not getting
> any error but there is no plot begin created. I will be grateful if someone
> could guide me about it. Here is the code:
>
> *********************************************************
> ;************************************************
> 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"
> ;************************************************
> ;Read data to plot and mask
> ;************************************************
> function mask_data_with_India_country(country_code,data)
> begin
>
> mask_start_time = get_cpu_time() ; For timing results
>
> ;---Convert rectilinear grid to 2D grid, but laid out as 1D array.
> dims = dimsizes(data)
> lat1d = ndtooned(conform_dims(dims,data&$data!0$,0))
> lon1d = ndtooned(conform_dims(dims,data&$data!1$,1))
>
> shapefile_name = "IND_adm0.shp"
> ;---Open shapefile and read lat/lon values.
> sfilename = "IND_adm0.shp"
> f = addfile(sfilename,"r")
> segments = f->segments
> geometry = f->geometry
> segsDims = dimsizes(segments)
> geomDims = dimsizes(geometry)
>
> ;---Read global attributes
> geom_segIndex = f at geom_segIndex
> geom_numSegs = f at geom_numSegs
> segs_xyzIndex = f at segs_xyzIndex
> segs_numPnts = f at segs_numPnts
> numFeatures = geomDims(0)
>
> lon = f->x
> lat = f->y
> nlatlon = dimsizes(lat)
>
> min_lat = min(lat)
> max_lat = max(lat)
> min_lon = min(lon)
> max_lon = max(lon)
>
> print("==================================================")
> print("Shapefile : " + sfilename)
> print("min_lat " + min_lat)
> print("max_lat " + max_lat)
> print("min_lon " + min_lon)
> print("max_lon " + max_lon)
>
> ii_latlon = ind(min_lat.le.lat1d.and.lat1d.le.max_lat.and.\
> min_lon.le.lon1d.and.lon1d.le.max_lon)
> nii = dimsizes(ii_latlon)
> print(nii + " values to check")
> print(numFeatures + " feature(s) to traverse with a maximum of " + \
> nlatlon + " points")
>
>
> ;---Create array to hold new data mask, and set all values to 0 initially.
> data_mask_1d = new(dimsizes(lat1d),integer)
> data_mask_1d = 0
>
> ;
> ; This is the loop that checks every point in lat1d/lon1d to see if it
> ; is inside or outside of the country. If it is inside, then data_mask_1d
> ; will be set to 1.
> ;
> ikeep = 0 ; Counter to see how many points were found inside the
> country
> do n=0,nii-1
> ii = ii_latlon(n)
> is_inside = False
> i = 0
> do while(.not.is_inside.and.i.lt.numFeatures)
> 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
> if(data_mask_1d(ii).ne.1.and.\
> gc_inout(lat1d(ii),lon1d(ii),\
> lat(startPT:endPT),lon(startPT:endPT))) then
> data_mask_1d(ii) = 1
> ikeep = ikeep+1
> is_inside = True
> continue
> end if
> end do
> i = i + 1
> end do
> end do
> print(ikeep + " values kept")
> print("==================================================")
>
> ; Create a 2D data array of same size as original data,
> ; but with appropriate values masked.
> ;
> data_mask = (where(onedtond(data_mask_1d,dims).eq.1,data,\
> data at _FillValue))
> copy_VarMeta(data,data_mask) ; Copy all metadata
>
> ;---Print timings
> mask_end_time = get_cpu_time()
> print("Elapsed time in CPU second for 'mask_data_with_India_country' = "
> + \
> (mask_end_time-mask_start_time))
>
> return(data_mask)
>
> ;----------------------------------------------------------------------
> ; Main code
> ;----------------------------------------------------------------------
> begin
> ;---Shapefile to use for masking. Using "Bolivia" here (BOL).
> India_code = "India"
> India_name = India_code + "_adm"
> India_shp = India_name + "/" + India_name + "0.shp"
>
> ;---Read lat/lon off shapefile
> s = addfile(India_shp,"r")
> slat = s->y
> slon = s->x
>
> ;---Read precipitation data to contour and mask
> f = addfile("b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECC.070001-
> 079912.nc","r")
> ts = f->PRECC(0,:,:)
> printVarSummary(ts)
>
> ;---Start the graphics
> wtype = "x11" ; "png"
> ; wtype at wkWidth = 2500 ; use for "png" or "x11"
> ; wtype at wkHeight = 2500
>
> wks = gsn_open_wks(wtype,"mask_India_"+India_code)
>
> res = True ; plot mods desired
>
> res at cnFillOn = True ; turn on color
> res at cnFillMode = "RasterFill"
> res at cnLinesOn = False ; turn off lines
> res at cnLineLabelsOn = False ; turn off labels
>
> res at mpMinLatF = min(ts&lat)
> res at mpMaxLatF = max(ts&lat)
> res at mpMinLonF = min(ts&lon)
> res at mpMaxLonF = max(ts&lon)
>
> ;---Create global plot of original data
> res at tiMainString = filename
> plot = gsn_csm_contour_map(wks,ts, res)
>
> ;---Mask "ts" against India shapefile outlines
> ts_mask = mask_data_with_India_country(India_code,ts)
> printVarSummary(ts_mask)
>
> ;---Set some additional resources for the second set of plots
>
> res at gsnDraw = False ; don't draw yet
> res at gsnFrame = False ; don't advance frame yet
>
> ;---Pick "nice" contour levels for both plots
> mnmxint = nice_mnmxintvl( min(ts_mask), max(ts_mask), 18, False)
> res at cnLevelSelectionMode = "ManualLevels"
> res at cnMinLevelValF = mnmxint(0)
> res at cnMaxLevelValF = mnmxint(1)
> res at cnLevelSpacingF = mnmxint(2)/2.
>
> res at mpMinLatF = min(slat)
> res at mpMaxLatF = max(slat)
> res at mpMinLonF = min(slon)
> res at mpMaxLonF = max(slon)
>
> res at gsnRightString = ""
> res at gsnLeftString = ""
>
> res at lbLabelBarOn = False
>
> ;---Create plot of original data
> res at tiMainString = "original data zoomed in"
> plot_orig = gsn_csm_contour_map(wks,ts, res)
>
> ;---Create plot of masked data
> res at tiMainString = "with country mask"
> plot_mask = gsn_csm_contour_map(wks,ts_mask, res)
>
> ;---Add shapefile outlines to both plots
> lnres = True
> lnres at gsLineColor = "NavyBlue"
>
> id_orig = gsn_add_shapefile_polylines(wks,plot_orig,India_shp,lnres)
> id_mask = gsn_add_shapefile_polylines(wks,plot_mask,India_shp,lnres)
>
> ;---Panel both plots on one page
> pres = True
> pres at txString = ts at long_name + " (" + ts at units + ")"
> pres at gsnMaximize = True
> pres at gsnPanelLabelBar = True
> gsn_panel(wks,(/plot_orig,plot_mask/),(/1,2/),pres)
> end
> end
>
>
> --
> ************************************************************
> *******************************************************
> "I slept and dreamt that life was joy. I awoke and saw that life was
> service. I acted and behold, service was joy." - Rabindranath Tagore
> ************************************************************
> ********************************************************
>
> Ipshita Majhi
> PhD Candidate
> University of Alaska , Fairbanks
> Atmospheric Science Department
> (907)978-4220 <(907)%20978-4220> ipmajhi at alaska.edu
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170308/47075c52/attachment.html
More information about the ncl-talk
mailing list