Как наложить shapely.Polygon на изображение в формате tiff ⇐ Python

[Anonymous]

У меня есть наклонное изображение здания, сохраненное в формате TIFF. Географическая привязка этого изображения осуществляется с помощью GCP (с использованием ESPG:32617). Кроме того, у меня есть shapely.Polygon, который представляет собой координаты одного и того же здания по широте и долготе (например, EPSG: 4326). Координаты и изображение из разных источников.
Я хочу отобразить и изображение, и многоугольник, используя одну и ту же систему отсчета, чтобы можно было проверить их совпадение.
Вот результаты выполнения gdalinfo в файле TIFF:

Драйвер: GTiff/GeoTIFF Файлы: oblique_N.tiff Размер 115, 317. Проекция GCP = PROJCRS["WGS 84 / зона UTM 17N", BASEGEOGCRS["WGS 84", ДАТУМ["Мировая геодезическая система 1984 года", ЭЛЛИПСОИД["WGS 84",6378137,298.257223563, LENGTHUNIT["метр",1]]], ПРИМЕМ["Гринвич",0, ANGLEUNIT["градус",0,0174532925199433]], Идентификатор["EPSG",4326]], CONVERSION["UTM зона 17N", МЕТОД["Поперечный Меркатор", Идентификатор["EPSG",9807]], ПАРАМЕТР["Широта естественного происхождения",0, ANGLEUNIT["градус",0,0174532925199433], Идентификатор["EPSG",8801]], ПАРАМЕТР["Долгота естественного происхождения",-81, ANGLEUNIT["градус",0,0174532925199433], Идентификатор["EPSG",8802]], ПАРАМЕТР["Масштабный коэффициент естественного происхождения",0,9996, SCALEUNIT["единство",1], Идентификатор["EPSG",8805]], ПАРАМЕТР["Ложное смещение",500000, LENGTHUNIT["метр",1], Идентификатор["EPSG",8806]], ПАРАМЕТР["Ложное северное смещение",0, LENGTHUNIT["метр",1], ID["EPSG",8807]]], CS[Декартово,2], ОСЬ["(E)",восток, ЗАКАЗ[1], LENGTHUNIT["метр",1]], ОСЬ["(N)",север, ЗАКАЗ[2], LENGTHUNIT["метр",1]], ИСПОЛЬЗОВАНИЕ[ ОБЛАСТЬ ПРИМЕНЕНИЯ["Инженерные изыскания, топографическое картографирование."], РАЙОН["Между 84° и 78° з.д., северное полушарие между экватором и 84° с.ш., на суше и на море. Багамские острова. Эквадор - к северу от экватора. Канада - Нунавут; Онтарио; Квебек. Каймановы острова. Колумбия. Коста-Рика. Куба Ямайка Никарагуа Панама. США (США)."], BBOX[0,-84,84,-78]], Идентификатор["EPSG",32617]] Сопоставление оси данных с осью CRS: 1,2 GCP[0]: идентификатор=1, информация= (98.3105276672193,0.120036521323868) -> (581069.63551352,2851093.99523912,271) GCP[1]: идентификатор=2, информация= (15.8348654576047,316.879963478676) -> (581056.130608335,2851042.88431689,271) GCP[2]: идентификатор=3, информация= (115304.815517974796) -> (581069.948644366,2851042.96910577,271) GCP[3]: идентификатор=4, информация= (0,15.4405225192958) -> (581055.817530873,2851093.91044916,271) Метаданные: AREA_OR_POINT=Площадь TIFFTAG_RESOLUTIONUNIT=1 (безразмерный) TIFFTAG_XRESOLUTION=1 TIFFTAG_YRESOLUTION=1 Метаданные структуры изображения: СЖАТИЕ=YCbCr JPEG МЕЖДУНИЕ=ПИКСЕЛЬ SOURCE_COLOR_SPACE=YCbCr Угловые координаты: Верхний левый (0,0, 0,0) Нижний левый (0,0, 317,0) Верхний правый ( 115.0, 0.0) Нижний правый ( 115.0, 317.0) Центр ( 57,5, 158,5) Блок 1 = 112x320 Тип = Байт, ColorInterp = Красный Нет данных Значение = 0 Блок 2 = 112x320 Тип = Байт, ColorInterp = Зеленый Нет данных Значение = 0 Блок 3 = 112x320 Тип = Байт, ColorInterp = Синий Нет данных Значение = 0 Вот многоугольник:

ПОЛИГОН ((-87.6388310614825 41.8948365045333, -87.6387464180662 41.8948376193115, -87.6386545 41.8948393, -87.6386589556933 41.8949865342181, -87.6388292759657 41.8949865342181, -87.6388310614825 41.8948365045333)) Вот моя попытка сделать это с помощью rasterio:

def show_tiff_oblique(имя файла, многоугольник, топор): src = rasterio.open(имя файла, 'r') tiff_crs = src.gcps[1] # Это EPSG:32617 polygon_crs = rasterio.crs.CRS.from_epsg(4326) src_transform = rasterio.transform.from_gcps(src.gcps[0]) # Вычислить преобразование вращения между заданным gcps и видом по умолчанию src_transform, ширина, высота = rasterio.warp.calculate_default_transform( tiff_crs, tiff_crs, src.width, src.height, gcps=src.gcps[0] ) # Проецируем полигон на CRS изображения projected_polygon = rasterio.warp.transform_geom(polygon_crs, tiff_crs,polygon.__geo_interface__) преобразователь = rasterio.transform.AffineTransformer(src_transform) очки = [] для x, y в projected_polygon['coordinates'][0]: Points.append(transformer.xy(x, y)) print((x, y), точки[-1]) reprojected_polygon = shapely.Polygon(точки) rasterio.plot.show(src, Transform=src_transform, ax=ax) ax.plot(*reprojected_polygon.exterior.xy, 'r') Полученный сюжет неверен, согласно нему здание и изображение находятся совершенно в разных местах, разница в миллионы.
Пример src_transform:

| 1,57, 0,00, 446931,48| | 0,00,-1,57, 4638425,53| | 0,00, 0,00, 1,00| Моя интуиция подсказывает мне, что я неправильно использую преобразования, в частности, мне кажется, что я неправильно получаю поворот изображения (это наклон, поэтому в большинстве случаев оно поворачивается). Однако мне не удалось найти пример того, как это получить.