[ncl-talk] Area average while using shapefile
Ipshita Majhi
ipmajhi at alaska.edu
Wed Mar 22 14:24:18 MDT 2017
Hi Mary,
I looked at the plot , it is creating a plot where it is masking value over
India, but I want value over India only with rest of the world masked out,I
am attaching the plot with this email.Could you suggest the required
changes. Thank you in advance.
On Fri, Mar 17, 2017 at 7:07 AM, Mary Haley <haley at ucar.edu> wrote:
> Hi Ipshita,
>
> The "conform" function doesn't do the masking. It simply propagates a 2D
> array to a 3D array.
>
> You need to take this 2D mask array, conform it to 3D, and then apply it
> to your original "ts" data in order to "mask" the locations in "ts" where
> the values are 1:
>
> ;---Create a 2D mask array using India shapefile outlines
> ts_mask_2d = mask_data_with_India_country(India_code,ts(0,:,:))
>
> ;---Conform this 2D array to a 3D array
> ts_mask_3d = conform(ts,ts_mask_2d,(/1,2/))
>
> ;---Apply the 3D mask to "ts"
> tm_3d = mask(ts,ts_mask_3d,1) ; Keep ts when ts_mask_3d equal 1
>
>
> The tm_3d array should now be the one that you can average.
>
> By the way, there are some things you can do to clean up your code a
> little.
>
> [1] You don't need to use a "do" loop to calculate the average. You can
> use "dim_avg_n" or "dim_avg_n_Wrap" (preserves metadata):
>
> prec_avg = dim_avg_n_Wrap(tm_3d,(/1,2/))
>
> [2] You generally don't ever need to use "new" unless you need to populate
> that array inside a do loop by indexing it, or you just need to create an
> empty array for some reason.
>
> So, you can replace this code:
>
> precp_mon_avg = new((/100,12/),"float")
>
> precp_mon_avg=reshape(prec_avg,(/100,12/))
>
> with just one line:
>
> precp_mon_avg = reshape(prec_avg,(/100,12/))
>
> --Mary
>
>
> On Thu, Mar 16, 2017 at 1:48 PM, Ipshita Majhi <ipmajhi at alaska.edu> wrote:
>
>> 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 <(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
>>
>>
>
--
*******************************************************************************************************************
"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/20170322/4f2f2f82/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mask.png
Type: image/png
Size: 134023 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170322/4f2f2f82/attachment.png
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PI_precc_sh_avg.ncl
Type: application/octet-stream
Size: 5426 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170322/4f2f2f82/attachment.obj
More information about the ncl-talk
mailing list