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

Question

我目前从GeoTiff文件中获得了一个171 x 171的图像（在其他情况下，我可能会有更大的图像），我的目标是将图像中的每个像素转换为纬度/经度对。我已经能够根据StackOverflow帖子将图像的角落转换为纬度/经度对：Obtain Latitude and Longitude from a GeoTIFF File，这篇帖子非常有帮助，因为我的原始坐标是在UTM区15中。

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

我看到的与此最相关的帖子是https://scriptndebug.wordpress.com/2014/11/24/latitudelongitude-of-each-pixel-using-python-and-gdal/。然而，该帖子建议基本上创建一个循环遍历每个像素并将其转换为纬度、经度。是否有更有效率的方法？

为了更好地说明我的实际用例，我的最终目标是有一堆卫星图像（例如，在此情况下，每个图像都是171 x 171）。我正试图创建一个建筑分割模型。现在，我正在尝试通过在每个图像上创建一个掩膜，将标记为1的像素与建筑物对应，否则标记为0，来产生带标签的数据点。首先，我使用Microsoft美国建筑物轮廓数据：https://github.com/microsoft/USBuildingFootprints，其中他们发布了由纬度、经度定义的多边形的GeoJSON文件，表示他们检测到的建筑物。我打算这样做的方式是：

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

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

- Vincent

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

- Vincent

---

这里导入geopandas和numba是不必要的。 - Dhyana