N-D indexing with defaults in NumPy

I don't know if there is anything in NumPy to do that directly, but you can always implement it yourself. This is not particularly smart or efficient, as it requires multiple advanced indexing operations, but does what you need:

import numpy as np

def get_with_default(a, indices, default=0):
    # Ensure inputs are arrays
    a = np.asarray(a)
    indices = tuple(np.broadcast_arrays(*indices))
    if len(indices) <= 0 or len(indices) > a.ndim:
        raise ValueError('invalid number of indices.')
    # Make mask of indices out of bounds
    mask = np.zeros(indices[0].shape, np.bool)
    for ind, s in zip(indices, a.shape):
        mask |= (ind < 0) | (ind >= s)
    # Only do masking if necessary
    n_mask = np.count_nonzero(mask)
    # Shortcut for the case where all is masked
    if n_mask == mask.size:
        return np.full_like(a, default)
    if n_mask > 0:
        # Ensure index arrays are contiguous so masking works right
        indices = tuple(map(np.ascontiguousarray, indices))
        for ind in indices:
            # Replace masked indices with zeros
            ind[mask] = 0
    # Get values
    res = a[indices]
    if n_mask > 0:
        # Replace values of masked indices with default value
        res[mask] = default
    return res

# Test
print(get_with_default(
    np.array([[1,2,3],[4,5,6]]),
    (np.array([0, 0, 1, 1, 2, 2]), np.array([1, 2, 2, 3, 3, 5])),
    13
))
# [ 2  3  6 13 13 13]

I arrived at this question because I was looking for exactly the same. I came up with the following function, which does what you ask for 2 dimension. It could likely be generalised to N dimensions.

def get_with_defaults(a, xx, yy, nodata):
   # get values from a, clipping the index values to valid ranges
   res = a[np.clip(yy, 0, a.shape[0] - 1), np.clip(xx, 0, a.shape[1] - 1)]
   # compute a mask for both x and y, where all invalid index values are set to true
   myy = np.ma.masked_outside(yy, 0, a.shape[0] - 1).mask
   mxx = np.ma.masked_outside(xx, 0, a.shape[1] - 1).mask
   # replace all values in res with NODATA, where either the x or y index are invalid
   np.choose(myy + mxx, [res, nodata], out=res)
   return res

xx and yy are the index array, a is indexed by (y,x).

This gives:

>>> a=np.zeros((3,2),dtype=int)
>>> get_with_defaults(a, (-1, 1000, 0, 1, 2), (0, -1, 0, 1, 2), -1)
array([-1, -1,  0,  0, -1])

As an alternative, the following implementation achieves the same and is more concise:

def get_with_default(a, xx, yy, nodata):
   # get values from a, clipping the index values to valid ranges
   res = a[np.clip(yy, 0, a.shape[0] - 1), np.clip(xx, 0, a.shape[1] - 1)]
   # replace all values in res with NODATA (gets broadcasted to the result array), where
   # either the x or y index are invalid
   res[(yy < 0) | (yy >= a.shape[0]) | (xx < 0) | (xx >= a.shape[1])] = nodata
   return res