Separating points on a map
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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment