Find Latitude/Longitude Coordinates of Every Pixel in a GeoTiff Image

Question

I currently have a 171 x 171 image from a GeoTiff file (although in other cases, I might have much bigger images). My goal is to take each pixel in the image and convert to latitude/longitude pair.

I am already able to convert the corners of the image to latitude/longitude pair based on this StackOverflow post: Obtain Latitude and Longitude from a GeoTIFF File. This post was helpful since my original coordinates were in UTM Zone 15.

However, I now want to convert all of the pixels of the image to latitude, longitude pair and store the results in a numpy array of the same dimension. So the output would be a numpy array that is 171 x 171 x 2 with each element of the numpy array being a tuple of the (longitude, latitude) pair.

The most relevant post I've seen on this is https://scriptndebug.wordpress.com/2014/11/24/latitudelongitude-of-each-pixel-using-python-and-gdal/. However, that post suggests to essentially create a for loop over each pixel and convert to latitude, longitude. Is there a way that is more efficient?

Just to give more context on my actual use case, my end goal is I have a bunch of satellite imagery (for example in this case, each image is 171 x 171). I am trying to create a building segmentation model. Right now, I am trying to produce labeled data points by creating a mask on each image that labels a pixel a 1 if it corresponds to a building, else 0. To start, I'm using the Microsoft US Building Footprint data: https://github.com/microsoft/USBuildingFootprints where they've released GeoJSON files of polygons (defined by latitude, longitude) of buildings they've detected. The way I'm thinking about doing this is:

1. Find the latitude, longitude of each pixel in my image. Thus, I will have 171 x 171 Points. Put this in a GeoSeries
2. Intersect the Points (in a GeoSeries) with the Microsoft US Building Footprints data (using GeoPandas intersects: https://geopandas.org/reference.html#geopandas.GeoSeries.intersects)
3. If the Point intersects with any Polygons in the Microsoft US Building Footprint data, then label 1, else 0.

Right now I'm on step (1), that is, efficiently find the latitude/longitude coordinate of each pixel in the image.

Answer 1

Unfortunately, I couldn't find a better solution (yet) than looping over all the pixels. Here's my solution so far:

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

Those were my helper functions. Because this would take a long time, in find_image_pixel_lat_lon_coord I saved the results into a dictionary image_coordinates_dict which I wrote to a pickle file to save results.

Then the way I would use this is:

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