Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
You can’t perform that action at this time.