[ncl-talk] Area average using shapefile

Ipshita Majhi ipmajhi at alaska.edu
Sat Mar 11 16:51:00 MST 2017

```Dear NCL,

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.

Is there a way to do the average for all the timestep?
Can I save the created mask and then use it like I do with landsea.nc and
use it for all diff precipitation dataset to get an average over India?

Here is the code that I am using to take an average

;************************************************
;************************************************
;************************************************
begin

mask_start_time = get_cpu_time()     ; For timing results

;---Convert rectilinear grid to 2D grid, but laid out as 1D array.
dims  = dimsizes(data)
lat1d = ndtooned(conform_dims(dims,data&\$data!0\$,0))
lon1d = ndtooned(conform_dims(dims,data&\$data!1\$,1))

;---Open shapefile and read lat/lon values.
segments = f->segments
geometry = f->geometry
segsDims = dimsizes(segments)
geomDims = dimsizes(geometry)

geom_segIndex = f at geom_segIndex
geom_numSegs  = f at geom_numSegs
segs_xyzIndex = f at segs_xyzIndex
segs_numPnts  = f at segs_numPnts
numFeatures   = geomDims(0)

lon  = f->x
lat  = f->y
nlatlon = dimsizes(lat)

min_lat = min(lat)
max_lat = max(lat)
min_lon = min(lon)
max_lon = max(lon)

print("==================================================")
print("Shapefile : "  + sfilename)
print("min_lat " + min_lat)
print("max_lat " + max_lat)
print("min_lon " + min_lon)
print("max_lon " + max_lon)

ii_latlon = ind(min_lat.le.lat1d.and.lat1d.le.max_lat.and.\
min_lon.le.lon1d.and.lon1d.le.max_lon)
nii = dimsizes(ii_latlon)
print(nii + " values to check")
print(numFeatures + " feature(s) to traverse with a maximum of " + \
nlatlon + " points")

;---Create array to hold new data mask, and set all values to 0 initially.

;
; This is the loop that checks every point in lat1d/lon1d to see if it
; is inside or outside of the country. If it is inside, then data_mask_1d
; will be set to 1.
;
ikeep = 0    ; Counter to see how many points were found inside the
country
do n=0,nii-1
ii = ii_latlon(n)
is_inside = False
i = 0
do while(.not.is_inside.and.i.lt.numFeatures)
startSegment = geometry(i, geom_segIndex)
numSegments  = geometry(i, geom_numSegs)
do seg=startSegment, startSegment+numSegments-1
startPT = segments(seg, segs_xyzIndex)
endPT   = startPT + segments(seg, segs_numPnts) - 1
gc_inout(lat1d(ii),lon1d(ii),\
lat(startPT:endPT),lon(startPT:endPT))) then
ikeep = ikeep+1
is_inside = True
continue
end if
end do
i = i + 1
end do
end do
print(ikeep + " values kept")
print("==================================================")

; Create a 2D data array of same size as original data,
; but with appropriate values masked.
;
data at _FillValue))

;---Print timings
print("Elapsed time in CPU second for 'mask_data_with_India_country' = "
+ \

end

;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
;---Shapefile to use for masking. Using "Bolivia" here (BOL).
India_code = "IND"
India_name = "IND" + "_adm0"

;---Read lat/lon off shapefile
slat      = s->y
slon      = s->x

;---Read precipitation data to contour and mask
filename = "b.e11.B1850C5CN.f09_g16.005.cam.h0.PRECSL.070001-079912.nc"
ts       = f->PRECSL
printVarSummary(ts)

;---Mask "ts" against India shapefile outlines
ts_mask at _FillValue = -2.1474*10^9

printVarSummary(data_avg)
print(data_avg)

end

--
*******************************************************************************************************************
"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
********************************************************************************************************************

Ipshita Majhi
PhD Candidate
University of Alaska , Fairbanks
Atmospheric Science Department
(907)978-4220 ipmajhi at alaska.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170311/80244882/attachment.html
```