<div dir="ltr">Dear NCL Users,<div><br></div><div>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:</div><div><br></div><div>*********************************************************</div><div><div>;************************************************</div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" </div><div>;************************************************</div><div>;Read data to plot and mask</div><div>;************************************************ </div><div>function mask_data_with_India_country(country_code,data) </div><div>begin</div><div><br></div><div>mask_start_time = get_cpu_time() ; For timing results</div><div><br></div><div>;---Convert rectilinear grid to 2D grid, but laid out as 1D array.</div><div> dims = dimsizes(data)</div><div> lat1d = ndtooned(conform_dims(dims,data&$data!0$,0))</div><div> lon1d = ndtooned(conform_dims(dims,data&$data!1$,1))</div><div><br></div><div> shapefile_name = "IND_adm0.shp"</div><div> ;---Open shapefile and read lat/lon values.</div><div> sfilename = "IND_adm0.shp"</div><div> f = addfile(sfilename,"r")</div><div> segments = f->segments</div><div> geometry = f->geometry</div><div> segsDims = dimsizes(segments)</div><div> geomDims = dimsizes(geometry)</div><div><br></div><div>;---Read global attributes </div><div> geom_segIndex = f@geom_segIndex</div><div> geom_numSegs = f@geom_numSegs</div><div> segs_xyzIndex = f@segs_xyzIndex</div><div> segs_numPnts = f@segs_numPnts</div><div> numFeatures = geomDims(0)</div><div> </div><div> lon = f->x</div><div> lat = f->y</div><div> nlatlon = dimsizes(lat)</div><div> </div><div> min_lat = min(lat)</div><div> max_lat = max(lat)</div><div> min_lon = min(lon)</div><div> max_lon = max(lon)</div><div><br></div><div> print("==================================================")</div><div> print("Shapefile : " + sfilename)</div><div> print("min_lat " + min_lat)</div><div> print("max_lat " + max_lat)</div><div> print("min_lon " + min_lon)</div><div> print("max_lon " + max_lon)</div><div> </div><div> ii_latlon = ind(min_lat.le.lat1d.and.lat1d.le.max_lat.and.\</div><div> min_lon.le.lon1d.and.lon1d.le.max_lon)</div><div> nii = dimsizes(ii_latlon)</div><div> print(nii + " values to check")</div><div> print(numFeatures + " feature(s) to traverse with a maximum of " + \</div><div> nlatlon + " points")</div><div> </div><div> </div><div> ;---Create array to hold new data mask, and set all values to 0 initially.</div><div> data_mask_1d = new(dimsizes(lat1d),integer)</div><div> data_mask_1d = 0</div><div><br></div><div>;</div><div>; This is the loop that checks every point in lat1d/lon1d to see if it</div><div>; is inside or outside of the country. If it is inside, then data_mask_1d</div><div>; will be set to 1.</div><div>;</div><div> ikeep = 0 ; Counter to see how many points were found inside the country</div><div> do n=0,nii-1</div><div> ii = ii_latlon(n)</div><div> is_inside = False</div><div> i = 0</div><div> do while(.not.is_inside.and.i.lt.numFeatures)</div><div> startSegment = geometry(i, geom_segIndex)</div><div> numSegments = geometry(i, geom_numSegs)</div><div> do seg=startSegment, startSegment+numSegments-1</div><div> startPT = segments(seg, segs_xyzIndex)</div><div> endPT = startPT + segments(seg, segs_numPnts) - 1</div><div> if(data_mask_1d(ii).ne.1.and.\</div><div> gc_inout(lat1d(ii),lon1d(ii),\</div><div> lat(startPT:endPT),lon(startPT:endPT))) then</div><div> data_mask_1d(ii) = 1</div><div> ikeep = ikeep+1</div><div> is_inside = True</div><div> continue</div><div> end if</div><div> end do</div><div> i = i + 1</div><div> end do</div><div> end do</div><div> print(ikeep + " values kept")</div><div> print("==================================================") </div><div><br></div><div>; Create a 2D data array of same size as original data,</div><div>; but with appropriate values masked.</div><div>;</div><div> data_mask = (where(onedtond(data_mask_1d,dims).eq.1,data,\</div><div> data@_FillValue))</div><div> copy_VarMeta(data,data_mask) ; Copy all metadata</div><div><br></div><div>;---Print timings</div><div> mask_end_time = get_cpu_time()</div><div> print("Elapsed time in CPU second for 'mask_data_with_India_country' = " + \</div><div> (mask_end_time-mask_start_time))</div><div><br></div><div> return(data_mask)</div><div> </div><div>;----------------------------------------------------------------------</div><div>; Main code</div><div>;----------------------------------------------------------------------</div><div>begin</div><div>;---Shapefile to use for masking. Using "Bolivia" here (BOL).</div><div>India_code = "India"</div><div>India_name = India_code + "_adm"</div><div>India_shp = India_name + "/" + India_name + "0.shp"</div><div><br></div><div>;---Read lat/lon off shapefile</div><div> s = addfile(India_shp,"r")</div><div> slat = s->y</div><div> slon = s->x</div><div><br></div><div>;---Read precipitation data to contour and mask</div><div> f = addfile("<a href="http://b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECC.070001-079912.nc">b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECC.070001-079912.nc</a>","r")</div><div> ts = f->PRECC(0,:,:)</div><div> printVarSummary(ts)</div><div><br></div><div>;---Start the graphics</div><div> wtype = "x11" ; "png"</div><div>; wtype@wkWidth = 2500 ; use for "png" or "x11"</div><div>; wtype@wkHeight = 2500</div><div><br></div><div> wks = gsn_open_wks(wtype,"mask_India_"+India_code)</div><div><br></div><div> res = True ; plot mods desired</div><div><br></div><div> res@cnFillOn = True ; turn on color</div><div> res@cnFillMode = "RasterFill"</div><div> res@cnLinesOn = False ; turn off lines</div><div> res@cnLineLabelsOn = False ; turn off labels</div><div><br></div><div> res@mpMinLatF = min(ts&lat)</div><div> res@mpMaxLatF = max(ts&lat)</div><div> res@mpMinLonF = min(ts&lon)</div><div> res@mpMaxLonF = max(ts&lon)</div><div><br></div><div>;---Create global plot of original data</div><div> res@tiMainString = filename</div><div> plot = gsn_csm_contour_map(wks,ts, res)</div><div><br></div><div>;---Mask "ts" against India shapefile outlines</div><div> ts_mask = mask_data_with_India_country(India_code,ts)</div><div> printVarSummary(ts_mask)</div><div><br></div><div>;---Set some additional resources for the second set of plots</div><div><br></div><div> res@gsnDraw = False ; don't draw yet</div><div> res@gsnFrame = False ; don't advance frame yet</div><div><br></div><div>;---Pick "nice" contour levels for both plots</div><div> mnmxint = nice_mnmxintvl( min(ts_mask), max(ts_mask), 18, False)</div><div> res@cnLevelSelectionMode = "ManualLevels"</div><div> res@cnMinLevelValF = mnmxint(0)</div><div> res@cnMaxLevelValF = mnmxint(1)</div><div> res@cnLevelSpacingF = mnmxint(2)/2.</div><div><br></div><div> res@mpMinLatF = min(slat)</div><div> res@mpMaxLatF = max(slat)</div><div> res@mpMinLonF = min(slon)</div><div> res@mpMaxLonF = max(slon)</div><div><br></div><div> res@gsnRightString = ""</div><div> res@gsnLeftString = ""</div><div><br></div><div> res@lbLabelBarOn = False</div><div><br></div><div>;---Create plot of original data</div><div> res@tiMainString = "original data zoomed in"</div><div> plot_orig = gsn_csm_contour_map(wks,ts, res)</div><div><br></div><div>;---Create plot of masked data</div><div> res@tiMainString = "with country mask"</div><div> plot_mask = gsn_csm_contour_map(wks,ts_mask, res)</div><div><br></div><div>;---Add shapefile outlines to both plots</div><div> lnres = True</div><div> lnres@gsLineColor = "NavyBlue"</div><div><br></div><div> id_orig = gsn_add_shapefile_polylines(wks,plot_orig,India_shp,lnres)</div><div> id_mask = gsn_add_shapefile_polylines(wks,plot_mask,India_shp,lnres)</div><div><br></div><div>;---Panel both plots on one page</div><div> pres = True</div><div> pres@txString = ts@long_name + " (" + ts@units + ")"</div><div> pres@gsnMaximize = True</div><div> pres@gsnPanelLabelBar = True</div><div> gsn_panel(wks,(/plot_orig,plot_mask/),(/1,2/),pres)</div><div>end</div><div>end</div><div><br></div><div><br></div>-- <br><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div><div><div>*******************************************************************************************************************</div><div><span style="font-size:small">"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</span><br></div><div><span style="font-size:small">********************************************************************************************************************</span></div><div><br></div><div>Ipshita Majhi<br></div>PhD Candidate<br></div>University of Alaska , Fairbanks<br></div>Atmospheric Science Department<br></div>(907)978-4220 <a href="mailto:ipmajhi@alaska.edu" target="_blank">ipmajhi@alaska.edu</a><br></div></div></div></div>
</div></div>