Reprojecting image with ground control points in geodetic CRS to equidistant cylindrical CRS on sphere

Question

I'm working with some satellite data from the Mars Reconnaissance Orbiter context camera (CTX). As a sanity check for some later processing that I would like to do, I'm trying to take the raster data I'm processing from the raw CTX experimental data records and georeference them against a ground truth product I found here. After fixing several issues, I've been able to get what's in the attached image. The transparent image on the top is mine and the background is the ground truth.

I believe the remaining error has something to do with the fact that the ground control points I am using on my image give geodetic lat and lon with respect to the Mars ellipsoid while the ground truth data product uses an equidistant cylindrical projection referenced to Mars sphere (specifically ESRI:103885 - Mars_2000_Equidistant_Cylindrical_sphere). This makes sense to me as the remaining error seems to be almost purely in the latitude direction and also appears to be greater closer to the +-45 degree line (at the bottom of the image here as these are in the southern hemisphere) which lines up with what I understand about how the error between the geodetic latitude and various auxiliary latitudes behaves.

However, I'm having trouble figuring out how I should deal with this discrepancy. I'm using rasterio to do the reprojection and I've reproduced the code below. Is there something I'm missing about how I should be preprocessing the latitude on my ground control points? Is there a way that I can get proj to handle this automatically?

# gcps is a list of ground control points with geodetic lat and lon in degrees        
gcps = self._ground_control_point_grid()

src_crs = CRS.from_proj4("+proj=longlat +a=3396190 +b=3376200 +no_defs")
dst_crs = CRS.from_proj4("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=3396190 +units=m +no_defs")

transform, width, height = calculate_default_transform(src_crs, dst_crs, self.width, self.height, gcps=gcps)

with rio.open(
    join(directory, fn),
    "w",
    driver="GTiff",
    height=height,
    width=width,
    count=1,
    dtype=self._view.dtype,
    crs=dst_crs,
    transform=transform,
    nodata=nodata
) as gtiff:
    reproject(self._view, rio.band(gtiff, 1), gcps=gcps, src_crs=src_crs, dst_crs=dst_crs, resampling=Resampling.nearest, num_threads=4, warp_mem_limit=1000, dst_nodata=nodata)