Check if any of the coordinates are within a certain distance using apply in Python pandas with two data frames
You need to drop your use of geopy. So far as I can see it does not support vectorised inverse geodesics. Instead use cartopy which inherently supports this.
The following example code uses random instead of real data, but with the same size as your real data, and completes in about one second. distances and close are both matrices whose rows correspond to vessels and columns correspond to ports. close is a matrix of booleans for vessel-port pairs that are proximate.
import numpy as np
from cartopy.geodesic import Geodesic
from numpy.random import default_rng

N_VESSELS = 4_800
N_PORTS = 137

rand = default_rng(seed=0)

vessels = rand.random((N_VESSELS, 2)) + (1, -51)
ports = rand.random((N_PORTS, 2)) + (1, -51)
source = vessels.repeat(N_PORTS, axis=0)
dest = np.tile(ports, (N_VESSELS, 1))

geo = Geodesic()
distances = geo.inverse(source, dest).base[:, 0].reshape((N_VESSELS, N_PORTS))
close = distances &lt; 3_000