Convert a GeoJSON polygon feature to a numpy array

    Args:
        feature (pygeoj.Feature): polygon feature to draw
        shape (tuple(int, int)): shape of 2D target numpy array to draw polygon in
        lat_idx (func): function converting a latitude to the (fractional) row index in the map
        lon_idx (func): function converting a longitude to the (fractional) column index in the map

    Returns:
        np.array: mask, background is zero, foreground is one

def geojson_polygon_to_mask(feature, shape, lat_idx, lon_idx):
    """Convert a GeoJSON polygon feature to a numpy array

    Args:
        feature (pygeoj.Feature): polygon feature to draw
        shape (tuple(int, int)): shape of 2D target numpy array to draw polygon in
        lat_idx (func): function converting a latitude to the (fractional) row index in the map
        lon_idx (func): function converting a longitude to the (fractional) column index in the map

    Returns:
        np.array: mask, background is zero, foreground is one
    """
    import matplotlib

    # specify 'agg' renderer, Mac renderer does not support what we want to do below
    matplotlib.use('agg')

    import matplotlib.pyplot as plt
    from matplotlib import patches
    import numpy as np

    # we can only do polygons right now
    if feature.geometry.type not in ('Polygon', 'MultiPolygon'):
        raise ValueError("Cannot handle feature of type " + feature.geometry.type)

    # fictional dpi - don't matter in the end
    dpi = 100

    # -- start documentation include: poly-setup
    # make a new figure with no frame, no axes, with the correct size, black background
    fig = plt.figure(frameon=False, dpi=dpi, )
    fig.set_size_inches(shape[1] / float(dpi), shape[0] / float(dpi))
    ax = plt.Axes(fig, [0., 0., 1., 1.])
    ax.set_axis_off()
    # noinspection PyTypeChecker
    ax.set_xlim([0, shape[1]])
    # noinspection PyTypeChecker
    ax.set_ylim([0, shape[0]])
    fig.add_axes(ax)
    # -- end documentation include: poly-setup

    # for normal polygons make coordinates iterable
    if feature.geometry.type == 'Polygon':
        coords = [feature.geometry.coordinates]
    else:
        coords = feature.geometry.coordinates

    for poly_coords in coords:
        # the polygon may contain multiple outlines; the first is
        # always the outer one, the others are 'holes'
        for i, outline in enumerate(poly_coords):
            # inside/outside fill value: figure background is white by
            # default, draw inverted polygon and invert again later
            value = 0. if i == 0 else 1.

            # convert lats/lons to row/column indices in the array
            outline = np.array(outline)
            xs = lon_idx(outline[:, 0])
            ys = lat_idx(outline[:, 1])

            # draw the polygon
            poly = patches.Polygon(list(zip(xs, ys)),
                                   facecolor=(value, value, value),
                                   edgecolor='none',
                                   antialiased=True)
            ax.add_patch(poly)

    # -- start documentation include: poly-extract
    # extract the figure to a numpy array,
    fig.canvas.draw()
    data = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
    # reshape to a proper numpy array, keep one channel only
    data = data.reshape(fig.canvas.get_width_height()[::-1] + (3,))[:, :, 0]
    # -- end documentation include: poly-extract

    # make sure we get the right shape back
    assert data.shape[0] == shape[0]
    assert data.shape[1] == shape[1]

    # convert from uints back to floats and invert to get black background
    data = 1. - data.astype(float) / 255.  # type: np.array

    # image is flipped horizontally w.r.t. map
    data = data[::-1, :]

    # done, clean up
    plt.close('all')

    return data