<div dir="ltr"><div>Dear NCL Users,</div><div><br></div><div>I am using shapefile to analyze data over India, thanks to Mary got the shape file working and then also got the data conformed. But when I am doing area average for each month the values are the same . I am not sure what is going wrong. I am pasting my code in this email.</div><div><br></div><div><br></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.shx"</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>end</div><div> </div><div>;----------------------------------------------------------------------</div><div>; Main code</div><div>;----------------------------------------------------------------------</div><div><br></div><div>;---Shapefile to use for masking. Using India (IND)</div><div>India_code = "IND"</div><div>India_name = "IND" + "_adm0"</div><div>India_shp = "IND_adm0.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> filename = "<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>"</div><div> f = addfile(filename,"r")</div><div> ts = f->PRECC(:,:,:)</div><div> </div><div><br></div><div>;---Mask "ts" against India shapefile outlines</div><div> ts_mask = mask_data_with_India_country(India_code,ts(0,:,:))</div><div> tm_3d=conform(ts,ts_mask,(/1,2/))</div><div><br></div><div> </div><div><br></div><div> tm_3d@_FillValue = ts@_FillValue</div><div> </div><div>; copy_VarAtts(ts,tm_3d) </div><div> copy_VarCoords(ts,tm_3d)</div><div><br></div><div> print(tm_3d(100,{25},{70}))</div><div> print(tm_3d(110,{25},{74}))</div><div><br></div><div> printVarSummary(tm_3d)</div><div> printMinMax(tm_3d, True)</div><div><br></div><div> prec_avg = new((/1200/),"float")</div><div><br></div><div> do ii =0,1199</div><div> </div><div> prec_avg(ii) = avg(tm_3d(ii,:,:))</div><div><br></div><div> end do</div><div><br></div><div> precp_mon_avg = new((/100,12/),"float")</div><div><br></div><div> precp_mon_avg=reshape(prec_avg,(/100,12/))</div><div><br></div><div><br></div><div> precp_mon = dim_avg_n(precp_mon_avg,0)</div><div><br></div><div> print(precp_mon)</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> </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>