def _landstat(landscape, updated_model, in_coords):
    """
    Compute the statistic for transforming coordinates onto an existing
    "landscape" of "mountains" representing source positions. Since the
    landscape is an array and therefore pixellated, the precision is limited.

    Parameters
    ----------
    landscape: nD array
        synthetic image representing locations of sources in reference plane
    updated_model: Model
        transformation (input -> reference) being investigated
    in_coords: nD array
        input coordinates

    Returns
    -------
    float:
        statistic representing quality of fit to be minimized
    """

    def _element_if_in_bounds(arr, index):
        try:
            return arr[index]
        except IndexError:
            return 0

    out_coords = updated_model(*in_coords)
    if len(in_coords) == 1:
        out_coords = (out_coords,)
    out_coords2 = tuple((coords - 0.5).astype(int) for coords in out_coords)
    result = sum(_element_if_in_bounds(landscape, coord[::-1]) for coord in zip(*out_coords2))
    ################################################################################
    # This stuff replaces the above 3 lines if speed doesn't hold up
    #    sum = np.sum(landscape[i] for i in out_coords if i>=0 and i<len(landscape))
    # elif len(in_coords) == 2:
    #    xt, yt = out_coords
    #    sum = np.sum(landscape[iy,ix] for ix,iy in zip((xt-0.5).astype(int),
    #                                                   (yt-0.5).astype(int))
    #                  if ix>=0 and iy>=0 and ix<landscape.shape[1]
    #                                     and iy<landscape.shape[0])
    ################################################################################
    return -result