How to overlay shapely.Polygon on tiff image

Question

I have an oblique image of a building, saved in a TIFF. That image is georeferenced using GCPs (using ESPG:32617). Additionally, I have a shapely.Polygon that's coordinates of the same building in latitude and longitude (so EPSG:4326). Coordinates and the image are from different sources.
I want to display both the image and polygon using the same reference frame, so I can verify that they align.
Here are results of gdalinfo on the TIFF file (link to file here):

Driver: GTiff/GeoTIFF
Files: oblique_N.tiff
Size is 115, 317
GCP Projection = 
PROJCRS["WGS 84 / UTM zone 17N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 17N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-81,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Between 84°W and 78°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Ecuador - north of equator. Canada - Nunavut; Ontario; Quebec. Cayman Islands. Colombia. Costa Rica. Cuba. Jamaica. Nicaragua. Panama. United States (USA)."],
        BBOX[0,-84,84,-78]],
    ID["EPSG",32617]]
Data axis to CRS axis mapping: 1,2
GCP[  0]: Id=1, Info=
          (98.3105276672193,0.120036521323868) -> (581069.63551352,2851093.99523912,271)
GCP[  1]: Id=2, Info=
          (15.8348654576047,316.879963478676) -> (581056.130608335,2851042.88431689,271)
GCP[  2]: Id=3, Info=
          (115,304.815517974796) -> (581069.948644366,2851042.96910577,271)
GCP[  3]: Id=4, Info=
          (0,15.4405225192958) -> (581055.817530873,2851093.91044916,271)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  COMPRESSION=YCbCr JPEG
  INTERLEAVE=PIXEL
  SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  317.0)
Upper Right (  115.0,    0.0)
Lower Right (  115.0,  317.0)
Center      (   57.5,  158.5)
Band 1 Block=112x320 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=112x320 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=112x320 Type=Byte, ColorInterp=Blue
  NoData Value=0

Here's the polygon:

POLYGON ((-87.6388310614825 41.8948365045333, -87.6387464180662 41.8948376193115, -87.6386545 41.8948393, -87.6386589556933 41.8949865342181, -87.6388292759657 41.8949865342181, -87.6388310614825 41.8948365045333))

Here's my attempt to do that using rasterio:

import rasterio
import rasterio.warp
import rasterio.plot
import matplotlib.pyplot as plt
import shapely


def show_tiff_oblique(filename, polygon, ax):
    def show_tiff_oblique(filename, polygon, ax):
    src = rasterio.open(filename, 'r')
    tiff_crs = src.gcps[1]
    polygon_crs = rasterio.crs.CRS.from_epsg(4326)
    src_transform = rasterio.transform.from_gcps(src.gcps[0])
    projected_polygon = rasterio.warp.transform_geom(polygon_crs, tiff_crs, polygon.__geo_interface__)
    
    points = []
    for x, y in projected_polygon['coordinates'][0]:
        points.append(~src_transform * (x, y))
    reprojected_polygon = shapely.Polygon(points)
    
    rasterio.plot.show(src, transform=src_transform, ax=ax)
    ax.plot(*reprojected_polygon.exterior.xy, 'r')

Resulting plot is wrong, according to it building and image are in totally different places, the difference is millions in magnitude. Resulting plot doesn't quite align with the building, sometimes quite wildy. I can't tell if it's due to the quality of data or if it's a matter of bad transformation somewhere along the way.
Sample ~src_transform:

| 0.74,-0.02,-259398.65|
| 0.49,-0.52, 2201773.68|
| 0.00, 0.00, 1.00|

My intuition tells me that I'm not using transforms correctly, in particular I think I'm not retrieving rotation of the image properly (it's an oblique so in most cases it's rotated). However, I wasn't able to find an example of how to retrieve that.

EDIT: Added link to the TIFF file.
EDIT2: Improved code based on my own findings.

Answer 1

This is not an answer to the question, it's just an analysis of the data.

First I verified the location of the .tiff in Qgis, as shown below. The area inside the red box is the .tiff image, and the lower layer is Google Satellite.The .tiff image is located in Miami.

Here we have a coordinate for example : x : 581063.887, y :2851069.988

Then I checked a coordinate in the Plygon, lat : 41.8948365045333 lon : -87.6388310614825, which corresponds this picture :

Which is located in Chicago.

The distance between the .tiff and the Polygon is nearly 1906Km.

So the coordinates of the polygon should be wrong and you need to reselect your polygon.