[ncl-talk] Area average while using shapefile
Ipshita Majhi
ipmajhi at alaska.edu
Fri Mar 17 11:50:49 MDT 2017
Thank you Mary for the solving the issue. I am thankful and grateful
Best Regards
Ipshita
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/20170317/c5787e8c/attachment.html
More information about the ncl-talk
mailing list