查找 GeoTiff 图像中每个像素的纬度/经度坐标

我目前有一个来自 GeoTiff 文件的 171 x 171 图像（尽管在其他情况下，我可能有更大的图像）。我的目标是获取图像中的每个像素并将其转换为纬度/经度对。

我已经能够根据此 StackOverflow 帖子将图像的角点转换为纬度/经度对：从 GeoTIFF 文件获取纬度和经度 https://stackoverflow.com/questions/2922532/obtain-latitude-and-longitude-from-a-geotiff-file。这篇文章很有帮助，因为我的原始坐标位于 UTM Zone 15。

但是，我现在想将图像的所有像素转换为纬度、经度对，并将结果存储在相同维度的 numpy 数组中。因此输出将是一个 171 x 171 x 2 的 numpy 数组，其中 numpy 数组的每个元素都是（经度，纬度）对的元组。

我见过的最相关的帖子是https://scriptndebug.wordpress.com/2014/11/24/latitudelongitude-of-each-pixel-using-python-and-gdal/ https://scriptndebug.wordpress.com/2014/11/24/latitudelongitude-of-each-pixel-using-python-and-gdal/。然而，该帖子建议本质上在每个像素上创建一个 for 循环并转换为纬度、经度。有没有更高效的方法呢？

为了提供有关我的实际用例的更多背景信息，我的最终目标是拥有一堆卫星图像（例如在本例中，每个图像都是 171 x 171）。我正在尝试创建一个建筑分割模型。现在，我尝试通过在每个图像上创建一个掩码来生成标记数据点，该掩码将像素标记为 1（如果对应于建筑物），否则标记为 0。 首先，我使用 Microsoft 美国建筑足迹数据：https://github.com/microsoft/USBuildingFootprints https://github.com/microsoft/USBuildingFootprints他们发布了检测到的建筑物的多边形（由纬度、经度定义）的 GeoJSON 文件。我考虑这样做的方式是：

1. 查找图像中每个像素的纬度和经度。因此，我将拥有 171 x 171 积分。将其放入 GeoSeries 中
2. 将点（在 GeoSeries 中）与 Microsoft 美国建筑足迹数据相交（使用 GeoPandas 相交：https://geopandas.org/reference.html#geopandas.GeoSeries.intersects https://geopandas.org/reference.html#geopandas.GeoSeries.intersects)
3. 如果该点与 Microsoft 美国建筑足迹数据中的任何多边形相交，则标记为 1，否则标记为 0。

现在我正在进行步骤（1），即有效地找到图像中每个像素的纬度/经度坐标。

不幸的是，我还找不到比循环所有像素更好的解决方案。到目前为止，这是我的解决方案：

import glob
import os
import pickle
import sys

import gdal
import geopandas as gpd
import matplotlib
import matplotlib.pyplot as plt
from numba import jit
import numpy as np
from osgeo import osr
import PIL
from PIL import Image, TiffImagePlugin
from shapely.geometry import Point, Polygon, box
import torch


def pixel2coord(img_path, x, y):
    """
    Returns latitude/longitude coordinates from pixel x, y coords

    Keyword Args:
      img_path: Text, path to tif image
      x: Pixel x coordinates. For example, if numpy array, this is the column index
      y: Pixel y coordinates. For example, if numpy array, this is the row index
    """
    # Open tif file
    ds = gdal.Open(img_path)

    old_cs = osr.SpatialReference()
    old_cs.ImportFromWkt(ds.GetProjectionRef())

    # create the new coordinate system
    # In this case, we'll use WGS 84
    # This is necessary becuase Planet Imagery is default in UTM (Zone 15). So we want to convert to latitude/longitude
    wgs84_wkt = """
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]]"""
    new_cs = osr.SpatialReference()
    new_cs.ImportFromWkt(wgs84_wkt)

    # create a transform object to convert between coordinate systems
    transform = osr.CoordinateTransformation(old_cs,new_cs) 
    
    gt = ds.GetGeoTransform()

    # GDAL affine transform parameters, According to gdal documentation xoff/yoff are image left corner, a/e are pixel wight/height and b/d is rotation and is zero if image is north up. 
    xoff, a, b, yoff, d, e = gt

    xp = a * x + b * y + xoff
    yp = d * x + e * y + yoff

    lat_lon = transform.TransformPoint(xp, yp) 

    xp = lat_lon[0]
    yp = lat_lon[1]
    
    return (xp, yp)


def find_img_coordinates(img_array, image_filename):
    img_coordinates = np.zeros((img_array.shape[0], img_array.shape[1], 2)).tolist()
    for row in range(0, img_array.shape[0]):
        for col in range(0, img_array.shape[1]): 
            img_coordinates[row][col] = Point(pixel2coord(img_path=image_filename, x=col, y=row))
    return img_coordinates


def find_image_pixel_lat_lon_coord(image_filenames, output_filename):
    """
    Find latitude, longitude coordinates for each pixel in the image

    Keyword Args:
      image_filenames: A list of paths to tif images
      output_filename: A string specifying the output filename of a pickle file to store results

    Returns image_coordinates_dict whose keys are filenames and values are an array of the same shape as the image with each element being the latitude/longitude coordinates.
    """
    image_coordinates_dict = {}
    for image_filename in image_filenames:
        print('Processing {}'.format(image_filename))
        img = Image.open(image_filename)
        img_array = np.array(img)
        img_coordinates = find_img_coordinates(img_array=img_array, image_filename=image_filename)
        image_coordinates_dict[image_filename] = img_coordinates
        with open(os.path.join(DATA_DIR, 'interim', output_filename + '.pkl'), 'wb') as f:
            pickle.dump(image_coordinates_dict, f)
    return image_coordinates_dict

这些是我的辅助功能。因为这需要很长时间find_image_pixel_lat_lon_coord我将结果保存到字典中image_coordinates_dict我将其写入 pickle 文件以保存结果。

那么我使用它的方式是：

# Create a list with all tif imagery
image_filenames = glob.glob(os.path.join(image_path_dir, '*.tif'))

image_coordinates_dict = find_image_pixel_lat_lon_coord(image_filenames, output_filename='image_coordinates')