[ncl-talk] Area average using shapefile

Mary Haley haley at ucar.edu
Sun Mar 12 13:43:47 MDT 2017


I think you just need to conform your ts_mask array to be the same size as
your data array with the multiple timesteps.

First, create the mask using a single timestep of "ts".  For example, if
"ts" is ntime x nlat x nlon:

  ts_mask = mask_data_with_India_country(India_code,ts(0,:,:))   ;
arbitrarily chose time index=0

It shouldn't matter which timestep you choose, as long as your lat/lon grid
is the same for every time step.

Now, conform ts_mask to be the same size as ts:

  ts_mask_3d = conform(ts,ts_mask,(/1,2/))

I think the rest of the code should work as expected.

--Mary



On Sat, Mar 11, 2017 at 4:51 PM, Ipshita Majhi <ipmajhi at alaska.edu> wrote:

> 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
>
> ;************************************************
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
> ;************************************************
> ;Read data to plot and mask
> ;************************************************
> function mask_data_with_India_country(country_code,data)
> 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))
>
>   shapefile_name = "IND_adm0.shx"
>  ;---Open shapefile and read lat/lon values.
>   sfilename = "IND_adm0.shp"
>   f = addfile(sfilename,"r")
>   segments = f->segments
>   geometry = f->geometry
>   segsDims = dimsizes(segments)
>   geomDims = dimsizes(geometry)
>
> ;---Read global attributes
>   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.
>   data_mask_1d = new(dimsizes(lat1d),integer)
>   data_mask_1d = 0
>
> ;
> ; 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
>          if(data_mask_1d(ii).ne.1.and.\
>             gc_inout(lat1d(ii),lon1d(ii),\
>                      lat(startPT:endPT),lon(startPT:endPT))) then
>            data_mask_1d(ii) = 1
>            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_mask = (where(onedtond(data_mask_1d,dims).eq.1,data,\
>               data at _FillValue))
>   copy_VarMeta(data,data_mask)      ; Copy all metadata
>
> ;---Print timings
>   mask_end_time = get_cpu_time()
>   print("Elapsed time in CPU second for 'mask_data_with_India_country' = "
> + \
>          (mask_end_time-mask_start_time))
>
>   return(data_mask)
> end
>
> ;----------------------------------------------------------------------
> ; Main code
> ;----------------------------------------------------------------------
> begin
> ;---Shapefile to use for masking. Using "Bolivia" here (BOL).
> India_code = "IND"
> India_name = "IND" + "_adm0"
> India_shp  = "IND_adm0.shp"
>
> ;---Read lat/lon off shapefile
>  s         = addfile(India_shp,"r")
>  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"
>   f        =  addfile(filename,"r")
>   ts       = f->PRECSL
>   printVarSummary(ts)
>
> ;---Mask "ts" against India shapefile outlines
>   ts_mask = mask_data_with_India_country(India_code,ts)
>   printVarSummary(ts_mask)
>  ts_mask at _FillValue = -2.1474*10^9
>  area_avg=wgt_areaave_Wrap(ts_mask,1,1,0)
>  data_avg  = wgt_areaave_Wrap((where(ts_mask.ge.0,ts_mask,ts at _FillValue)),1,1,0)
>
>  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 <(907)%20978-4220> ipmajhi at alaska.edu
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170312/07c22160/attachment.html 


More information about the ncl-talk mailing list