<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">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>;************************************************<br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" <br>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl" <br>;************************************************<br>;Read data to plot and mask<br>;************************************************ <br>function mask_data_with_India_country(country_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,data&$data!0$,0))<br> lon1d = ndtooned(conform_dims(dims,data&$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("==================================================")<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.le.max_lat.and.\<br> min_lon.le.lon1d.and.lon1d.le.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(.not.is_inside.and.i.lt.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: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("==================================================") <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,dims).eq.1,data,\<br> data@_FillValue))<br> copy_VarMeta(data,data_mask) ; 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))<br><br> return(data_mask)<br>end<br> <br>;----------------------------------------------------------------------<br>; Main code<br>;----------------------------------------------------------------------<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">b.e11.B1850C5CN.f09_g16.005.cam.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(India_code,ts)<br> printVarSummary(ts_mask)<br> ts_mask@_FillValue = -2.1474*10^9<br> area_avg=wgt_areaave_Wrap(ts_mask,1,1,0)<br> data_avg = wgt_areaave_Wrap((where(ts_mask.ge.0,ts_mask,ts@_FillValue)),1,1,0) <br> printVarSummary(data_avg)<br>print(data_avg)<br><br><br><br><br>end<br clear="all"><div><div><div><div><div><div><br>-- <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 target="_blank" href="mailto:ipmajhi@alaska.edu">ipmajhi@alaska.edu</a><br></div></div></div></div>
</div></div></div></div></div></div></div>