查找GeoTiff图像中每个像素的纬度/经度坐标 [英] Find Latitude/Longitude Coordinates of Every Pixel in a GeoTiff Image

问题描述

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

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.

我已经能够基于以下StackOverflow帖子将图像的角点转换为纬度/经度对:

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.

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

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?

仅在实际使用情况下提供更多背景信息，我的最终目标是获得一堆卫星图像(例如，在这种情况下，每个图像为171 x 171).我正在尝试创建建筑物分割模型.现在，我试图通过在每个图像上创建一个遮罩来产生标记的数据点，该遮罩将像素标记为1(如果它对应于建筑物)，否则标记为0.首先，我使用的是Microsoft US Building Footprint数据: https://github.com/microsoft/USBuildingFootprints ，他们已发布多边形的GeoJSON文件(由纬度定义)，经度).我正在考虑这样做的方式是:

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. 查找图像中每个像素的纬度和经度.因此，我将得到171 x 171点.将此放入GeoSeries
2. 将点(在GeoSeries中)与Microsoft美国建筑足迹数据相交(使用GeoPandas相交:现在我进入步骤(1)，也就是说，有效地找到图像中每个像素的纬度/经度坐标.

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

   推荐答案

   不幸的是，我还没有找到比遍历所有像素更好的解决方案.到目前为止，这是我的解决方案:

   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
   

   这些是我的助手功能.因为这会花费很长时间，所以在 find_image_pixel_pixel_lat_lon_coord 中，我将结果保存到字典 image_coordinates_dict 中，并写入了一个pickle文件中以保存结果.

   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.

   那我使用它的方式是:

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

   这篇关于查找GeoTiff图像中每个像素的纬度/经度坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！