[ncl-talk] Area average while using shapefile
Ipshita Majhi
ipmajhi at alaska.edu
Fri Mar 24 14:00:55 MDT 2017
Thank you Mary for all the suggestions and support.
Best Regards
Ipshita
On Fri, Mar 24, 2017 at 7:20 AM, Mary Haley <haley at ucar.edu> wrote:
> Hi Ipshita,
>
> Sorry for misspelling your name earlier.
>
> I was able to run your script using a previous data file you sent me, and
> I realize that your "mask" array is not set up the way I thought.
>
> Your original mask command was very close, but remember what I said about
> what the function does. It actually *keeps* (protects) the values of your
> first array at locations where your second array equals the third array.
>
> So, in your case, you were effectively keeping all the values of precc at
> locations where precc_mask3d was equal to -9999, which is the opposite of
> what you want.
>
> A better way to think of this is "at locations where precc_mask_3d is not
> missing, keep the values of precc":
>
> precc_3d = mask(precc,.not.ismissing(precc_mask_3d),True)
> You could also use the where function:
>
> precc_3d = where(.not.ismissing(precc_mask_3d),precc,precc at _FillValue)
>
> --Mary
>
>
>
>
> On Thu, Mar 23, 2017 at 12:36 PM, Ipshita Majhi <ipmajhi at alaska.edu>
> wrote:
>
>> Hi Mary,
>>
>> I have tried that in the past as well. I am not sure what I am missing
>> but , it is not able to plot if I put 1 instead of -9999
>> It gives the following error:-
>> warning:ContourPlotIntialize: no valid values in scalar field;
>> ContourPlot not possible:[errno=1101]
>> I am not sure how to resolve it. I will be grateful for any suggestions
>>
>> On Thu, Mar 23, 2017 at 6:09 AM, Mary Haley <haley at ucar.edu> wrote:
>>
>>> 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.c
>>>>>> am.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
>>>>
>>>
>>>
>>
>>
>> --
>> ************************************************************
>> *******************************************************
>> "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
>>
>
>
--
*******************************************************************************************************************
"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/20170324/6881363d/attachment.html
More information about the ncl-talk
mailing list