[ncl-talk] Regarding masking with shape file in NCL

Pawan kumar chaubey acppawan at gmail.com
Thu May 5 23:01:28 MDT 2022


Dear All,
I am Facing an error while masking the data with the shapefile.
Please guide me and resolve my problem.


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

  load "./shapefile_utils.ncl"
  load "./mask_data_with_gadm_country.ncl"

  diri   = "./"
  fili   = "IMD_2000_2010_daily_rf_0.25_ncl_monsum.nc"



  f      = addfile(diri+fili, "r")

  prc    = flt2dble(f->rf)
  pmsg   = prc at _FillValue      ; convenience

  printVarSummary(prc)
  printMinMax(prc,0)

  runlen = (/ 12 /)
  nrun   = dimsizes(runlen)
;*********************************
; plot parameters
;*********************************
  dimprc = dimsizes(prc)
  ntim   = dimprc(0)
  nlat   = dimprc(1)
  mlon   = dimprc(2)


  yyyymm = f->time
  yyyymm = yyyymm/100
  year    = yyyymm/100
  yrStrt  = toint(year(0))
  yrLast  = toint(year(ntim-1))
  nyear   = yrLast-yrStrt+1
  yyyymm  = yyyymm_time(yrStrt, yrLast, "integer")
  yrfrac  =  (/yyyymm_to_yyyyfrac(yyyymm, 0.0)/)
  prc&time = (/yyyymm/)


;=================================
wks_type = "png"
   wks_type at wkWidth = 2500      ; for resolution
   wks_type at wkHeight = 2500
  wks          = gsn_open_wks ("png","spi_mask") ; send graphics to PNG file
  cmap = read_colormap_file("BlueWhiteOrangeRed")  ; read color data

  res                      = True
  res at gsnDraw              = False              ; don't draw
  res at gsnFrame             = False              ; don't advance frame
  res at gsnAddCyclic         = False
  res at cnFillOn             = True               ; color fill
  res at cnFillPalette        = cmap(::-1,:)       ; set and reverse color map
  res at cnFillMode           = "RasterFill"       ; Raster Mode
  res at cnLinesOn            =  False             ; Turn off contour lines
  res at cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res at cnMinLevelValF       = -3.0               ; set min contour level
  res at cnMaxLevelValF       =  3.0               ; set max contour level
  res at cnLevelSpacingF      =  0.5               ; set contour spacing
  res at lbLabelBarOn         = False              ; turn off individual cb's
  res at txFont               = "courier-bold"     ; bold text
  res at mpFillOn             = False              ; backgroung map

  resP                      = True                   ; panel resources
 ;resP at gsnPanelMainString  = "SPI: 1980-2019 (IMD)"  ; add center string
  resP at gsnPanelLabelBar     = True                   ; add common colorbar

  res at mpMinLatF             =  7.5         ; India limits
  res at mpMaxLatF             =  37.5
  res at mpMinLonF             =  68
  res at mpMaxLonF             =  97.5



   diri   = "./"
 shp  = "BGD_adm1.shp"

 s      = addfile(diri+shp, "r")

 slat      = s->y
 slon      = s->x


 opt=True
 opt at spi_type = 3


  plot = new (1, "graphic")

  do nr=0,nrun-1
    spi    = dim_spi_n(prc, runlen(nr), False, 0)

   ; spi at long_name = "SPI"
   ; spi at units = "run="+runlen(nr)
    copy_VarCoords(prc, spi)
spi_mask   = shapefile_mask_data(spi,shp,opt)
end do



   opt         = True
 ; opt at minlat  = min(spi&latitude)
 ; opt at maxlat  = max(spi&latitude)
 ; opt at minlon  = min(spi&longitude)
 ; opt at maxlon  = max(spi&longitude)
 ; spi_mask    = shapefile_mask_data(spi,shp,opt)



  res at mpMinLatF       := min(slat)
  res at mpMaxLatF       := max(slat)
  res at mpMinLonF       := min(slon)
  res at mpMaxLonF       := max(slon)


  ;  res at gsnCenterString  = "1901-1910"

  ;  plot(0) = gsn_csm_contour_map(wks,spi_mask({201012},:,:), res)

  ;  gsn_panel(wks,plot,(/1,1/),resP)          ; now draw as one plot
;end do


;--------------------------------
; Create plot of masked data
;--------------------------------
  plot = gsn_csm_contour_map(wks,spi_mask({195012},:,:), res)
;end do
;---Add shapefile outlines to plots

  lnres             = True
  lnres at gsLineColor = "NavyBlue"
  id_mask           = gsn_add_shapefile_polylines(wks,plot_mask,shp,lnres)


draw(plot_mask)
frame(wks)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


-- 
|| with regards ||

*Pawan*
Research Scholar
Banaras Hindu University, Varanasi , India
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20220506/b3fbeaa8/attachment.html>


More information about the ncl-talk mailing list