[ncl-talk] Area average while using shapefile
Mary Haley
haley at ucar.edu
Thu Mar 23 08:09:19 MDT 2017
Ipsita,
I think the problem is with the "mask" function.
I always have a hard time remembering how this function works, so I
remember this:
The mask function protects the first argument at locations where the
second argument is equal to the third argument.
The script has this line:
precc_3d=mask(precc,precc_mask_3d,-9999)
Which effectively is going to keep the values of "precc" at locations where
precc_mask_3d is equal to -9999.
I think what it should be is:
precc_3d=mask(precc,precc_mask_3d,1)
I'm assuming that "precc_mask_3d" is set to the value 1 at locations where
the data is inside India.
--Mary
On Wed, Mar 22, 2017 at 2:24 PM, Ipshita Majhi <ipmajhi at alaska.edu> wrote:
> 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 <(907)%20978-4220> ipmajhi at alaska.edu
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170323/cdf268d6/attachment.html
More information about the ncl-talk
mailing list