[ncl-talk] Mask using shapefile
Ipshita Majhi
ipmajhi at alaska.edu
Wed Mar 8 12:59:14 MST 2017
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 ipmajhi at alaska.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170308/e2d787f4/attachment.html
More information about the ncl-talk
mailing list