<div dir="ltr"><div class="gmail_default" style="font-size:small">Hi,</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Can you provide the b.e11 NetCDF file? You can email me offline about this if you want.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I do recognize the code you included because it's code that I wrote, but it's hard to debug without the data file.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">Also, did you look at the "shapefile_mask_data" function on the shapefiles example page?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">In particular, you could look at example shapefile_11.ncl:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default"><a href="http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex11">http://www.ncl.ucar.edu/Applications/shapefiles.shtml#ex11</a><br></div><div class="gmail_default"><br></div><div class="gmail_default">You can first try it with minimal options to see what happens:</div><div class="gmail_default"><br></div><div class="gmail_default"><div class="gmail_default"><font face="monospace, monospace"> opt = True</font></div><div class="gmail_default"><font face="monospace, monospace"> opt@debug = True<br></font></div><div class="gmail_default"><font face="monospace, monospace"> ts_mask = shapefile_mask_data(ts,India_shp,opt)<br> printVarSummary(ts_mask)</font></div><div class="gmail_default"><font face="monospace, monospace"> printMinMax(ts_mask,0)</font></div><div class="gmail_default"><br></div><div class="gmail_default">I can't remember, but you might need to call this too:</div><div class="gmail_default"><br></div><div class="gmail_default"><div style="font-size:12.8px"><font face="monospace, monospace"> copy_VarMeta(ts,ts_mask) ; Copy all metadata</font></div><div><br></div></div></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">--Mary</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Mar 8, 2017 at 12:59 PM, Ipshita Majhi <span dir="ltr"><<a href="mailto:ipmajhi@alaska.edu" target="_blank">ipmajhi@alaska.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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>******************************<wbr>***************************</div><div><div>;*****************************<wbr>*******************</div><div>load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_code.ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_csm.ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/contributed.<wbr>ncl" </div><div>load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/shea_util.ncl" </div><div>;*****************************<wbr>*******************</div><div>;Read data to plot and mask</div><div>;*****************************<wbr>******************* </div><div>function mask_data_with_India_country(<wbr>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,<wbr>data&$data!0$,0))</div><div> lon1d = ndtooned(conform_dims(dims,<wbr>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("=======================<wbr>===========================")</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.<wbr>lat1d.le.max_lat.and.\</div><div> min_lon.le.lon1d.and.lon1d.le.<wbr>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(.<a href="http://not.is_inside.and.i.lt">not.is_inside.and.i.lt</a>.<wbr>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.<wbr>\</div><div> gc_inout(lat1d(ii),lon1d(ii),\</div><div> lat(startPT:endPT),lon(<wbr>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("=======================<wbr>===========================") </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,<wbr>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_<wbr>time))</div><div><br></div><div> return(data_mask)</div><div> </div><div>;-----------------------------<wbr>------------------------------<wbr>-----------</div><div>; Main code</div><div>;-----------------------------<wbr>------------------------------<wbr>-----------</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" target="_blank">b.e11.B1850C5CN.f09_<wbr>g16.005.cam.h0.PRECC.070001-<wbr>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_<wbr>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(<wbr>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_<wbr>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(<wbr>wks,plot_orig,India_shp,lnres)</div><div> id_mask = gsn_add_shapefile_polylines(<wbr>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,<wbr>plot_mask/),(/1,2/),pres)</div><div>end</div><div>end</div><span class="HOEnZb"><font color="#888888"><div><br></div><div><br></div>-- <br><div class="m_-7272105917876497920gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div><div><div>******************************<wbr>******************************<wbr>******************************<wbr>*************************</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">******************************<wbr>******************************<wbr>******************************<wbr>**************************</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><a href="tel:(907)%20978-4220" value="+19079784220" target="_blank">(907)978-4220</a> <a href="mailto:ipmajhi@alaska.edu" target="_blank">ipmajhi@alaska.edu</a><br></div></div></div></div>
</font></span></div></div>
<br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>