Why are my "loadings" offset in my PCA biplot? (Python, Matplotlib)

I know there are packages that will do biplots for me but I want to understand how these loadings are plotted because as of now it seems like wizardry.

My understanding is the following:

1. Get the loadings
2. Calculate the relative contribution of each feature in each PC
3. Use the contributions as magnitudes, then use the directionality from the original loadings
4. Plot those vectors

I'm comparing my results against the PCA package (https://github.com/erdogant/pca) as I trust these results more than my intuition. I'm trying to figure where my logic is flawed and what I'm doing wrong with my plotting. The original source code starts plotting the biplot here.

The iris example looks accurate but my small sample coral gene expression climate change dataset doesn't look like the biplots are correct at all.

I just can't seem to figure out what I'm doing wrong. Any help would be great.

import matplotlib.pyplot as plt
from matplotlib.patches import ConnectionPatch
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler

# Test version to understand how loadings are plotted
def plot_biplot(
    data:pd.DataFrame,
    arrow_scale=1,
    n_feat=4,
    ):

    # Build PCA
    model = PCA(n_components=min(data.shape))
    projection = pd.DataFrame(model.fit_transform(data), index=data.index)
    projection.columns = projection.columns.map(lambda j: "PC.{}".format(j+1))
    projection = projection.iloc[:,:2]

    # Get loading subset
    loadings = pd.DataFrame(model.components_.T, index=data.columns)
    loadings.columns = loadings.columns.map(lambda j: "PC.{}".format(j+1))
    loadings = loadings.iloc[:,:2]
    features = loadings.iloc[:,0].abs().sort_values(ascending=False).iloc[:n_feat].index
    loadings = loadings.loc[features]

    # Get contributions of each feature to each PC
    contributions = loadings.abs()
    contributions = contributions/contributions.sum(axis=0)
    contributions = contributions.iloc[:,:2]

    # Get projections
    x = projection.iloc[:,0].values
    y = projection.iloc[:,1].values

    mean_x = np.mean(x)
    max_x = np.abs(x).max()
    mean_y = np.mean(y)
    max_y = np.abs(y).max()

    # Keywords
    _arrow_kws = dict(color="red", linewidth=1, alpha=0.618) 
    _arrowhead_kws = dict(arrowstyle="-|>", mutation_scale=20, facecolor="red", alpha=0.618)
    _text_kws = dict(color="red", ha='center', va='center')

    # Plot arrows and text
    with plt.style.context("seaborn-white"):
        fig, ax = plt.subplots(figsize=(5,5))
        ax.scatter(x,y,s=50,c="black", linewidth=1, edgecolor="white")
        ax.set_xlabel("PC.1({:.3f})".format(model.explained_variance_ratio_[0]), fontsize=15)
        ax.set_ylabel("PC.2({:.3f})".format(model.explained_variance_ratio_[1]), fontsize=15)

        for i in range(0, loadings.shape[0]):
            feature = loadings.index[i]
            # Set PC1 vs PC2 direction. Note that these are not neccarily the best loading.
            xarrow = contributions.loc[feature].iloc[0] * np.sign(loadings.loc[feature].iloc[0]) * max_x * arrow_scale
            yarrow = contributions.loc[feature].iloc[1] * np.sign(loadings.loc[feature].iloc[1]) * max_y * arrow_scale

            # Plot arrow
            xyA = (mean_x, mean_y)
            xyB = (xarrow, yarrow)
            con = ConnectionPatch(xyA, xyB, coordsA="data", coordsB="data", **_arrowhead_kws)
            ax.plot([xyA[0], xyB[0]], [xyA[1], xyB[1]], **_arrow_kws)
            ax.add_artist(con)

            # Plot feature label
            ax.text(xarrow*1.11, yarrow*1.11, feature,  **_text_kws)

        return ax
    
# Gold standard
# !pip install pca
# https://github.com/erdogant/pca
from pca import pca

def pca_wrapper(data, n_feat=4):
    # Initialize to reduce the data up to the number of componentes that explains 95% of the variance.
    model = pca(n_components=0.95)

    # Or reduce the data towards 2 PCs
    model = pca(n_components=2)

    # Fit transform
    results = model.fit_transform(data)

    # Make biplot with the number of features
    with plt.style.context("seaborn-white"):
        fig, ax = model.biplot(y=np.asarray(["_"]*data.shape[0]), n_feat=n_feat, legend=False, label=False, figsize=(5,5))
    return ax

# Iris data
iris = load_iris()
X_iris = pd.DataFrame(
    data=iris.data, 
    columns=map(lambda j:j[:-5].replace(" ","_"), iris.feature_names),
    index=map(lambda i:"iris_{}".format(i), range(150))
)

# Test version
plot_biplot(X_iris, arrow_scale=2)

# Gold standard
pca_wrapper(X_iris)

# Coral dataset
X_corals = pd.read_csv("https://pastebin.com/raw/iVinU2p3", sep="\t", index_col=0)


# Test version
plot_biplot(X_corals, arrow_scale=2)

# Gold standard
pca_wrapper(X_corals)