[ncl-talk] Area average while using shapefile

[Ipshita Majhi]

Dear NCL Users,

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.


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

;---Shapefile to use for masking. Using India (IND)
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.PRECC.070001-079912.nc"
  f        =  addfile(filename,"r")
  ts       = f->PRECC(:,:,:)


;---Mask "ts" against India shapefile outlines
  ts_mask = mask_data_with_India_country(India_code,ts(0,:,:))
  tm_3d=conform(ts,ts_mask,(/1,2/))



  tm_3d at _FillValue = ts at _FillValue

;  copy_VarAtts(ts,tm_3d)
  copy_VarCoords(ts,tm_3d)

 print(tm_3d(100,{25},{70}))
 print(tm_3d(110,{25},{74}))

 printVarSummary(tm_3d)
 printMinMax(tm_3d, True)

 prec_avg = new((/1200/),"float")

 do ii =0,1199

 prec_avg(ii) = avg(tm_3d(ii,:,:))

 end do

 precp_mon_avg = new((/100,12/),"float")

 precp_mon_avg=reshape(prec_avg,(/100,12/))


 precp_mon = dim_avg_n(precp_mon_avg,0)

 print(precp_mon)








-- 
*******************************************************************************************************************
"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/20170316/64093a20/attachment.html