"""
Module contain functions to visualize/plot the carthographic data.
"""
import logging
from typing_extensions import Literal
import numpy as np
import pandas as pd
from scipy import interpolate
from cmocean import cm
import matplotlib.pyplot as plt
[docs]def plot_nautical_chart(
coords: pd.DataFrame, ax: plt.Axes = None, mark_depths: bool = False
):
"""
Takes a pandas dataframe with x,y,z coordinates and creates a nautical chart based
on a linear interpolation of the specified coordinates. Uses the topo colormap from cmocean
Parameters
----------
coords:
Coordinates of measured depths should contain columns named "x", "y", "z".
ax:
the matplotlib axis to plot on. If no axis is specified, a new figure will be created.
mark_depths:
if true plots the coordinates on which the interpolation is based.
Returns
-------
f, ax
the figure and axis object
"""
if ax is None:
f, ax = plt.subplots(figsize=(20, 15))
# ax.pcolormesh(X,Y,Z)
interpolated_dataset = interpolate_heatmap(
coords["x"], coords["y"], coords["z"], n=400
)
levels = np.arange(-500, 201, 10)
ax.contour(
interpolated_dataset[0],
interpolated_dataset[1],
interpolated_dataset[2],
levels=levels,
colors="black",
linestyles="solid",
linewidths=0.5,
)
from matplotlib import colors
divnorm = colors.TwoSlopeNorm(vmin=-500.0, vcenter=0.0, vmax=200)
mappable = ax.contourf(
interpolated_dataset[0],
interpolated_dataset[1],
interpolated_dataset[2],
cmap=cm.topo,
levels=levels,
norm=divnorm,
# vmin=-500,
# vcenter=0,
# vmax=500,
extend="min",
)
# mappable.set_clim(-500, 500)
cbar = f.colorbar(mappable, ax=ax)
# cbar.
# Scatter plot of measured datapoints
if mark_depths:
ax.plot(
coords["x"],
coords["y"],
marker=".",
linewidth=0.5,
linestyle="",
markerfacecolor="white",
markeredgecolor="red",
markersize=3,
c="white",
alpha=1,
)
ax.set_aspect("equal")
ax.grid()
ax.set_xlabel("Longtitude (m)")
ax.set_ylabel("Lattitude (m)")
cbar.set_label("Elevation (m)")
return f, ax
def areas(ip):
p = ip.tri.points[ip.tri.vertices]
q = p[:, :-1, :] - p[:, -1, None, :]
areas = abs(q[:, 0, 0] * q[:, 1, 1] - q[:, 0, 1] * q[:, 1, 0]) / 2
return areas
def scale(points, xy_mean, xy_scale):
points = np.asarray(points, dtype=float)
return (points - xy_mean) / xy_scale
def unscale(points, xy_mean, xy_scale):
points = np.asarray(points, dtype=float)
return points * xy_scale + xy_mean
[docs]def interpolate_heatmap(
x,
y,
z,
n: int = None,
interp_method: Literal["linear", "nearest"] = "linear",
):
"""
The output of this method can directly be used for plt.imshow(z_grid, extent=extent, aspect='auto')
where the extent is determined by the min and max of the x_grid and y_grid.
The output can also be used as input for ax.pcolormesh(x, y, Z,**kw)
Parameters
----------
x : :class:`numpy.ndarray`
x data points
y : :class:`numpy.ndarray`
y data points
z : :class:`numpy.ndarray`
z data points
n
number of points for each dimension on the interpolated grid if set to None will auto determine amount of
points needed
interp_method
determines what interpolation method is used.
Returns
-------
x_grid : :class:`numpy.ndarray`
N*1 array of x-values of the interpolated grid
y_grid : :class:`numpy.ndarray`
N*1 array of x-values of the interpolated grid
z_grid : :class:`numpy.ndarray`
N*N array of z-values that form a grid.
"""
points = list(zip(x, y))
lbrt = np.min(points, axis=0), np.max(points, axis=0)
lbrt = lbrt[0][0], lbrt[0][1], lbrt[1][0], lbrt[1][1]
xy_mean = np.mean([lbrt[0], lbrt[2]]), np.mean([lbrt[1], lbrt[3]])
xy_scale = np.ptp([lbrt[0], lbrt[2]]), np.ptp([lbrt[1], lbrt[3]])
# interpolation needs to happen on a rescaled grid, this is somewhat akin
# to an assumption in the interpolation that the scale of the experiment
# is chosen sensibly. N.B. even if interp_method == "nearest" the linear
# interpolation is used to determine the amount of grid points.
ip = interpolate.LinearNDInterpolator(
scale(points, xy_mean=xy_mean, xy_scale=xy_scale), z
)
if n is None:
# Calculate how many grid points are needed.
# factor from A=√3/4 * a² (equilateral triangle)
# N.B. a factor 4 was added as there were to few points for uniform grid otherwise.
all_areas = areas(ip)
area_min = all_areas[all_areas > 0.0].min()
n = int(0.658 / np.sqrt(area_min)) * 4
n = max(n, 10)
if n > 500:
logging.debug("n: {} larger than 500".format(n))
n = 500
x_lin = y_lin = np.linspace(-0.5, 0.5, n)
if interp_method == "linear":
z_grid = ip(x_lin[:, None], y_lin[None, :]).squeeze()
elif interp_method == "nearest":
ip = interpolate.NearestNDInterpolator(
scale(points, xy_mean=xy_mean, xy_scale=xy_scale), z
)
z_grid = ip(x_lin[:, None], y_lin[None, :]).squeeze()
else:
raise ValueError("interp_method must be one of {'linear', 'nearest', 'deg'}")
# x and y grid points need to be rescaled from the linearly chosen points
points_grid = unscale(list(zip(x_lin, y_lin)), xy_mean=xy_mean, xy_scale=xy_scale)
x_grid = points_grid[:, 0]
y_grid = points_grid[:, 1]
return x_grid, y_grid, z_grid.T