Separating points on a map

separate_points.py

def lat_lon_hex_mesh(bounds, d=3):
    """Creates a hexagonal lat/lon mesh within bounds that has a radial separation of d"""

    lone, lonw, lats, latn = bounds

    # heigt of equilatral triangle from radial distance sqrt(r^2 - (r/2)^2)
    h = np.sqrt(0.75) * d
    w = d / 2

    lat_vals = np.arange(lats, latn, h)
    lon_vals = np.arange(lone, lonw, w)

    data = np.tile([[0, 1], [1, 0]],
                   (int(np.ceil(len(lat_vals) / 2)),
                    int(np.ceil(len(lon_vals) / 2))))
    if len(lat_vals) % 2 != 0 or len(lon_vals) % 2 != 0:
        data = data[0:len(lat_vals), 0:len(lon_vals)]

    try:
        grid = pd.DataFrame(data, index=lat_vals, columns=lon_vals)
    except:
        import ipdb
        ipdb.set_trace()

    return grid


def align_points_to_grid(df, bounds, d, grid_func=lat_lon_hex_mesh):

    grid = grid_func(bounds, d)

    new_df = df[['lat', 'lon']].copy(deep=True)
    new_df['key'] = [base64.b64encode(s.encode()).decode()[::-1] for s in df.index]
    new_df.sort_values('key')  # ensure psuedo-random mapping, to minimize pushing

    for site, point in new_df[['lat', 'lon']].iterrows():
        if point.lon < bounds[0] or bounds[1] < point.lon or \
                point.lat < bounds[2] or bounds[3] < point.lat:
            # skip sites outside of bounds
            continue
        left = grid.columns.searchsorted(point.lon - 4 * d)
        right = grid.columns.searchsorted(point.lon + 4 * d, side='right')
        bottom = grid.index.searchsorted(point.lat - 4 * d)
        top = grid.index.searchsorted(point.lat + 4 * d, side='right')
        close_grid_points = grid.iloc[bottom:top, left:right]

        best_dist = np.inf
        grid_loc = [np.nan, np.nan]
        for lat, col in close_grid_points.iterrows():
            for lon, val in col.iteritems():
                dist = vincenty(tuple(point), (lat, lon)).km
                if val == 1 and dist < best_dist:
                    best_dist = dist
                    grid_loc = [lat, lon]

        if np.isnan(grid_loc[0]):
            import ipdb
            ipdb.set_trace()
            assert False, 'Failed to find a close empty grid-point - try a finer grid, or a larger radius?'

        grid.loc[grid_loc[0], grid_loc[1]] = 2
        new_df.loc[site, ['lat', 'lon']] = grid_loc

    return new_df


def p_basic_map(df, bounds, ax, d=0):

    df = df.copy(deep=True)

    lonw, lone, lats, latn = bounds
    m = Basemap(ax=ax, llcrnrlon=lonw, llcrnrlat=lats, urcrnrlon=lone, urcrnrlat=latn)
    m.drawmapboundary()  # fill_color='aqua')
    m.fillcontinents()  # color='coral', lake_color='aqua')
    # m.drawcoastlines()

    # draw parallels and meridians.
    # label parallels on right and top
    # meridians on bottom and left
    parallels = np.arange(-90., 91, 30.)
    # labels = [left, right, top, bottom]
    m.drawparallels(parallels, labels=[False, True, True, False],
                    dashes=[100, 0.0001], color='grey')
    meridians = np.arange(0., 361., 30.)
    m.drawmeridians(meridians, labels=[True, False, False, True],
                    dashes=[100, 0.0001], color='grey')

    if d > 0:
        new_lat_lon = align_points_to_grid(df, bounds, d=d)

        for i in range(len(new_lat_lon)):
            m.plot([df.iloc[i]['lon'], new_lat_lon.iloc[i]['lon']],
                   [df.iloc[i]['lat'], new_lat_lon.iloc[i]['lat']],
                   c='0.5')

        points = m.scatter(x=new_lat_lon['lon'], y=new_lat_lon['lat'], c=df['c'],
                           cmap=get_segmented_colourmap(df['c']),
                           latlon=True, zorder=10,
                           edgecolors='1', lw=0.01)
    else:
        points = m.scatter(x=df['lon'], y=df['lat'], c=df['c'],
                           cmap=get_segmented_colourmap(df['c']),
                           latlon=True, zorder=10,
                           edgecolors='1', lw=0.01)

    return m, points