Gist for the following Stack Overflow question: https://stackoverflow.com/questions/54545758/create-equal-aspect-square-plot-with-multiple-axes-when-data-limits-are-differ

squares.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as grid

from mpl_toolkits.axes_grid1 import make_axes_locatable, axes_size


def joint_plot(x, y, ax=None):
    """
    Create a square joint plot of x and y.
    """
    if ax is None:
        ax = plt.gca()

    divider = make_axes_locatable(ax)
    xhax = divider.append_axes("top", size=1, pad=0.1, sharex=ax)
    yhax = divider.append_axes("right", size=1, pad=0.1, sharey=ax)

    ax.scatter(x, y)
    xhax.hist(x)
    yhax.hist(y, orientation="horizontal")

    x0,x1 = ax.get_xlim()
    y0,y1 = ax.get_ylim()
    ax.set_aspect(abs(x1-x0)/abs(y1-y0))

    plt.sca(ax)
    return ax, xhax, yhax


def color_plot(x, y, colors, ax=None):
    if ax is None:
        ax = plt.gca()

    divider = make_axes_locatable(ax)
    cbax = divider.append_axes("right", size="5%", pad=0.1)

    sc = ax.scatter(x, y, marker='o', c=colors, cmap='RdBu')
    plt.colorbar(sc, cax=cbax)

    ax.set_aspect("equal")

    plt.sca(ax)
    return ax, cbax


if __name__ == "__main__":
    _, axes = plt.subplots(nrows=2, ncols=2, figsize=(9,6))

    # Plot 1
    x = np.random.normal(100, 17, 120)
    y = np.random.normal(0.5, 0.1, 120)
    joint_plot(x, y, axes[0,0])

    # Plot 2
    x = np.random.normal(100, 17, 120)
    y = np.random.normal(100, 17, 120)
    c = np.random.normal(100, 17, 120)
    color_plot(x, y, c, axes[0,1])

    # Plot 3
    x = np.random.normal(100, 17, 120)
    y = np.random.normal(0.5, 0.1, 120)
    c = np.random.uniform(0.0, 1.0, 120)
    color_plot(x, y, c, axes[1,0])

    # Plot 4
    x = np.random.normal(0.5, 0.1, 120)
    y = np.random.normal(0.5, 0.1, 120)
    joint_plot(x, y, axes[1,1])

    plt.tight_layout()
    plt.show()