I have the following sample code (pretty simple algorithm that uses a rolling filter window) and am wondering what the best way is of speeding it up.  I tried rewriting it in Cython by pre-declaring the variables but that didn’t buy me a lot of time.  Then I rewrote it in Fortran (and compiled it with f2py) and now it’s lightning fast. 

But I would still like to know if I could rewrite it in pure python/numpy/scipy

or in Cython and get a similar speedup.

 

Here is the raw Python code: 

 

def mixed_coastline_slow(nsidc, radius, count, mask=None):

 

    nsidc_copy = numpy.copy(nsidc)

 

    if (mask is None):

        idx_coastline = numpy.where(nsidc_copy == NSIDC_COASTLINE_MIXED)

    else:

        idx_coastline = numpy.where(mask & (nsidc_copy == NSIDC_COASTLINE_MIXED))

       

    for (irow0, icol0) in zip(idx_coastline[0], idx_coastline[1]):

 

        rows = ( max(irow0-radius, 0), min(irow0+radius+1, nsidc_copy.shape[0]) )

        cols = ( max(icol0-radius, 0), min(icol0+radius+1, nsidc_copy.shape[1]) )

        window = nsidc[rows[0]:rows[1], cols[0]:cols[1]]

 

        npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()

        nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \

                                True, False).sum()

 

        if (100.0*nsnowice/npoints >= count):

             nsidc_copy[irow0, icol0] = MISR_SEAICE_THRESHOLD

 

    return nsidc_copy

 

and here is my attempt at Cython-izing it:

 

import numpy

cimport numpy as cnumpy

cimport cython

 

cdef int NSIDC_SIZE  = 721

cdef int NSIDC_NO_SNOW = 0

cdef int NSIDC_ALL_SNOW = 100

cdef int NSIDC_FRESHSNOW = 103

cdef int NSIDC_PERMSNOW  = 101

cdef int NSIDC_SEAICE_LOW  = 1

cdef int NSIDC_SEAICE_HIGH = 100

cdef int NSIDC_COASTLINE_MIXED = 252

cdef int NSIDC_SUSPECT_ICE = 253

 

cdef int MISR_SEAICE_THRESHOLD = 6

 

def mixed_coastline(cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc, int radius, int count):

 

     cdef int irow, icol, irow1, irow2, icol1, icol2, npoints, nsnowice

     cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc2 \

        = numpy.empty(shape=(NSIDC_SIZE, NSIDC_SIZE), dtype=numpy.uint8)

     cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] window \

        = numpy.empty(shape=(2*radius+1, 2*radius+1), dtype=numpy.uint8)

 

     nsidc2 = numpy.copy(nsidc)

 

     idx_coastline = numpy.where(nsidc2 == NSIDC_COASTLINE_MIXED)

       

     for (irow, icol) in zip(idx_coastline[0], idx_coastline[1]):

 

          irow1 = max(irow-radius, 0)

          irow2 = min(irow+radius+1, NSIDC_SIZE)

          icol1 = max(icol-radius, 0)

          icol2 = min(icol+radius+1, NSIDC_SIZE)

          window = nsidc[irow1:irow2, icol1:icol2]

 

          npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()

          nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \

                                  True, False).sum()

 

          if (100.0*nsnowice/npoints >= count):

               nsidc2[irow, icol] = MISR_SEAICE_THRESHOLD

 

     return nsidc2

 

Thanks in advance for any advice!

 

Catherine

 

_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@python.org
https://mail.python.org/mailman/listinfo/numpy-discussion