[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