[Numpy-discussion] subsampling array without loops

[Catherine Moroney]

This almost works.  Is there a way to do some masking on tiles, for
instance taking the maximum height of each 2x2 square that is an
odd number?   I've tried playing around with masking and where, but
they don't return an array of the original size and shape of "tiles"
below.

Catherine

> Perhaps you could take the relevant slices of the arrays with  
> appropriate striding:
>
> a = numpy.arange(128*512).reshape((128, 512))
> top_left = a[::2, ::2]
> top_right = a[::2, 1::2]
> bottom_left = a[1::2, ::2]
> bottom_right = a[1::2, 1::2]
>
> tiles = numpy.array([top_left, top_right, bottom_left, bottom_right])
> maxes = tiles.max(axis=0)
>
> Similarly if you want overlapping tiles, you could leave out the  
> final :2 in the slice specifications above.
>
> Zach
>
>
>
>> I'm looking for a way to acccomplish the following task without lots
>> of loops involved, which are really slowing down my code.
>>
>> I have a 128x512 array which I want to break down into 2x2 squares.
>> Then, for each 2x2 square I want to do some simple calculations
>> such as finding the maximum value, which I then store in a 64x256
>> array.
>>
>> Here is the actual code involved.  It's only slightly more complex
>> than what I described above, since it also involves doing a bit of
>> masking on the 2x2 sub-arrays.
>>
>> Any hints are much appreciated!  An inordinate amount of time is
>> being spent in this function and another one like it.
>>
>> Catherine
>>
>>     def calc_sdcm_at_rlra(self,iblock):
>>
>>         npixl = self.nlineht/self.nlinerl
>>         npixs = self.nsmpht/self.nsmprl
>>
>>         sdcm_out = numpy.array([Constants.CLOUDMASK_NR] \
>>                                *self.nlinerl*self.nsmprl,'int8') \
>>                         .reshape(self.nlinerl,self.nsmprl)
>>
>>         for ilinerl in range(0,self.nlinerl):
>>             for ismprl in range(0,self.nsmprl):
>>
>>                 height = self.data[iblock].height 
>> [ilinerl*2:ilinerl*2+2, \
>>                                                  ismprl*2:ismprl*2+2]
>>                 sdcm   = self.data[iblock].sdcm[ilinerl*2:ilinerl*2 
>> +2, \
>>                                                 ismprl*2:ismprl*2+2]
>>                 source = self.data[iblock].heightsrc 
>> [ilinerl*2:ilinerl*2+2, \
>>                                                       
>> ismprl*2:ismprl*2+2]
>>
>>                 mask1 = (source == Constants.HEIGHT_STEREO)
>>                 mask2 = ( (source == Constants.HEIGHT_SURFACE) | \
>>                           (source == Constants.HEIGHT_DEFAULT) )
>>
>>                 if (mask1.any()):
>>                     loc = height[mask1].argmax()
>>                     sdcm_out[ilinerl,ismprl] = sdcm[mask1].ravel() 
>> [loc]
>>                 elif (mask2.any()):
>>                     loc = height[mask2].argmax()
>>                     sdcm_out[ilinerl,ismprl] = sdcm[mask2].ravel() 
>> [loc]
>>
>>         return sdcm_out
>