如何在tiff图像上叠加shapely.Polygon

Question

我有一个建筑物的倾斜图像，保存在 TIFF 格式中。该图像使用 GCP 进行地理配准（使用 ESPG:32617）。此外，我有一个

shapely.Polygon

，它是同一建筑物的纬度和经度坐标（因此 EPSG：4326）。坐标和图像来自不同的来源。
我想使用相同的参考系显示图像和多边形，这样我就可以验证它们是否对齐。
以下是 TIFF 文件上

gdalinfo

的结果（链接到文件此处）：

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

这是多边形：

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))

这是我尝试使用光栅来做到这一点：

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')

结果图是错误的，根据它，建筑物和图像完全不同的地方，差异是数百万级的。 由此产生的情节与建筑不太相符，有时甚至相当疯狂。我无法判断这是否是由于数据质量造成的，或者是否是过程中某个地方转换不当的问题。
样本

~src_transform

：

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

我的直觉告诉我，我没有正确使用变换，特别是我认为我没有正确检索图像的旋转（它是倾斜的，所以在大多数情况下它是旋转的）。但是，我无法找到如何检索它的示例。

编辑：添加了指向 TIFF 文件的链接。
EDIT2：根据我自己的发现改进了代码。

Answer 1

这不是问题的答案，只是数据分析。

首先我验证了Qgis中

.tiff

的位置，如下图。红框内的区域是.tiff图像，下层是Google Satellite。

.tiff

图像位于Miami。

这里我们有一个坐标例如：

x : 581063.887, y :2851069.988

然后我检查了

Plygon

、

lat : 41.8948365045333 lon : -87.6388310614825

中的一个坐标，对应的是这张图：

位于芝加哥。

.tiff

与

Polygon

之间的距离近1906Km。

所以多边形的坐标应该是错误的，你需要重新选择你的多边形。