<div dir="ltr">Thank you Mary, it worked and I am grateful your suggestion<div><br></div><div>Best Regards</div><div>Ipshita</div></div><div class="gmail_extra"><br><div class="gmail_quote">On Sun, Mar 12, 2017 at 11:43 AM, Mary Haley <span dir="ltr"><<a href="mailto:haley@ucar.edu" target="_blank">haley@ucar.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"><div class="gmail_default" style="font-size:small">I think you just need to conform your ts_mask array to be the same size as your data array with the multiple timesteps.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">First, create the mask using a single timestep of "ts". For example, if "ts" is ntime x nlat x nlon:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"> ts_mask = mask_data_with_India_country(</span><span style="font-size:12.8px">I<wbr>ndia_code,ts(0,:,:)) ; arbitrarily chose time index=0</span><br style="font-size:12.8px"></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px">It shouldn't matter which timestep you choose, as long as your lat/lon grid is the same for every time step.</span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px">Now, conform ts_mask to be the same size as ts:</span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"> ts_mask_3d = conform(ts,ts_mask,(/1,2/)) </span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px">I think the rest of the code should work as expected.</span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px">--Mary</span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div><div class="gmail_default" style="font-size:small"><span style="font-size:12.8px"><br></span></div></div><div class="gmail_extra"><br><div class="gmail_quote"><div><div class="h5">On Sat, Mar 11, 2017 at 4:51 PM, Ipshita Majhi <span dir="ltr"><<a href="mailto:ipmajhi@alaska.edu" target="_blank">ipmajhi@alaska.edu</a>></span> wrote:<br></div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div class="h5"><div dir="ltr"><div><div><div><div>Dear NCL,<br><br></div>I got the shapefile plot running thanks to Mary. I wanted to take an average of precipitaiton over India, if I run the code for one time step I am able to take an average but if I put try to do it for all the timesteps it gives me error.<br><br></div>Is there a way to do the average for all the timestep?<br></div>Can I save the created mask and then use it like I do with <a href="http://landsea.nc" target="_blank">landsea.nc</a> and use it for all diff precipitation dataset to get an average over India?<br><br></div>Here is the code that I am using to take an average<br><br>;*****************************<wbr>*******************<br>load "$NCARG_ROOT/lib/ncarg/nclscri<wbr>pts/csm/gsn_code.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscri<wbr>pts/csm/gsn_csm.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscri<wbr>pts/csm/contributed.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscri<wbr>pts/csm/shea_util.ncl" <br>;*****************************<wbr>*******************<br>;Read data to plot and mask<br>;*****************************<wbr>******************* <br>function mask_data_with_India_country(c<wbr>ountry_code,data) <br>begin<br><br>mask_start_time = get_cpu_time() ; For timing results<br><br>;---Convert rectilinear grid to 2D grid, but laid out as 1D array.<br> dims = dimsizes(data)<br> lat1d = ndtooned(conform_dims(dims,dat<wbr>a&$data!0$,0))<br> lon1d = ndtooned(conform_dims(dims,dat<wbr>a&$data!1$,1))<br><br> shapefile_name = "IND_adm0.shx"<br> ;---Open shapefile and read lat/lon values.<br> sfilename = "IND_adm0.shp"<br> f = addfile(sfilename,"r")<br> segments = f->segments<br> geometry = f->geometry<br> segsDims = dimsizes(segments)<br> geomDims = dimsizes(geometry)<br><br>;---Read global attributes <br> geom_segIndex = f@geom_segIndex<br> geom_numSegs = f@geom_numSegs<br> segs_xyzIndex = f@segs_xyzIndex<br> segs_numPnts = f@segs_numPnts<br> numFeatures = geomDims(0)<br> <br> lon = f->x<br> lat = f->y<br> nlatlon = dimsizes(lat)<br> <br> min_lat = min(lat)<br> max_lat = max(lat)<br> min_lon = min(lon)<br> max_lon = max(lon)<br><br> print("=======================<wbr>===========================")<br> print("Shapefile : " + sfilename)<br> print("min_lat " + min_lat)<br> print("max_lat " + max_lat)<br> print("min_lon " + min_lon)<br> print("max_lon " + max_lon)<br> <br> ii_latlon = ind(min_lat.le.lat1d.and.lat1d<wbr>.le.max_lat.and.\<br> min_lon.le.lon1d.and.lon1d.le.<wbr>max_lon)<br> nii = dimsizes(ii_latlon)<br> print(nii + " values to check")<br> print(numFeatures + " feature(s) to traverse with a maximum of " + \<br> nlatlon + " points")<br> <br> <br> ;---Create array to hold new data mask, and set all values to 0 initially.<br> data_mask_1d = new(dimsizes(lat1d),integer)<br> data_mask_1d = 0<br><br>;<br>; This is the loop that checks every point in lat1d/lon1d to see if it<br>; is inside or outside of the country. If it is inside, then data_mask_1d<br>; will be set to 1.<br>;<br> ikeep = 0 ; Counter to see how many points were found inside the country<br> do n=0,nii-1<br> ii = ii_latlon(n)<br> is_inside = False<br> i = 0<br> do while(.<a href="http://not.is_inside.and.i.lt" target="_blank">not.is_inside.and.i.lt</a>.<wbr>numFeatures)<br> startSegment = geometry(i, geom_segIndex)<br> numSegments = geometry(i, geom_numSegs)<br> do seg=startSegment, startSegment+numSegments-1<br> startPT = segments(seg, segs_xyzIndex)<br> endPT = startPT + segments(seg, segs_numPnts) - 1<br> if(data_mask_1d(ii).ne.1.and.\<br> gc_inout(lat1d(ii),lon1d(ii),\<br> lat(startPT:endPT),lon(startPT<wbr>:endPT))) then<br> data_mask_1d(ii) = 1<br> ikeep = ikeep+1<br> is_inside = True<br> continue<br> end if<br> end do<br> i = i + 1<br> end do<br> end do<br> print(ikeep + " values kept")<br> print("=======================<wbr>===========================") <br><br>; Create a 2D data array of same size as original data,<br>; but with appropriate values masked.<br>;<br> data_mask = (where(onedtond(data_mask_1d,d<wbr>ims).eq.1,data,\<br> data@_FillValue))<br> copy_VarMeta(data,data_mask) <wbr> ; Copy all metadata<br><br>;---Print timings<br> mask_end_time = get_cpu_time()<br> print("Elapsed time in CPU second for 'mask_data_with_India_country' = " + \<br> (mask_end_time-mask_start_time<wbr>))<br><br> return(data_mask)<br>end<br> <br>;-----------------------------<wbr>------------------------------<wbr>-----------<br>; Main code<br>;-----------------------------<wbr>------------------------------<wbr>-----------<br>begin<br>;---Shapefile to use for masking. Using "Bolivia" here (BOL).<br>India_code = "IND"<br>India_name = "IND" + "_adm0"<br>India_shp = "IND_adm0.shp"<br><br>;---Read lat/lon off shapefile<br> s = addfile(India_shp,"r")<br> slat = s->y<br> slon = s->x<br><br>;---Read precipitation data to contour and mask<br> filename = "<a href="http://b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECSL.070001-079912.nc" target="_blank">b.e11.B1850C5CN.f09_g16.005.c<wbr>am.h0.PRECSL.070001-079912.nc</a>"<br> f = addfile(filename,"r")<br> ts = f->PRECSL<br> printVarSummary(ts)<br><br>;---Mask "ts" against India shapefile outlines<br> ts_mask = mask_data_with_India_country(I<wbr>ndia_code,ts)<br> printVarSummary(ts_mask)<br> ts_mask@_FillValue = -2.1474*10^9<br> area_avg=wgt_areaave_Wrap(ts_<wbr>mask,1,1,0)<br> data_avg = wgt_areaave_Wrap((where(ts_mas<wbr>k.ge.0,ts_mask,ts@_FillValue))<wbr>,1,1,0) <br> printVarSummary(data_avg)<br>print(data_avg)<br><br><br><br><br>end<span class="m_-7689327612168260257HOEnZb"><font color="#888888"><br clear="all"><div><div><div><div><div><div><br>-- <br><div class="m_-7689327612168260257m_-634011611513420806gmail_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>
</div></div></div></div></div></div></font></span></div>
<br></div></div>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">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/mailma<wbr>n/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="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>