使用 GDAL 将经度、纬度转换为像素值

Question

我有一个卫星 GeoTIFF 图像和一个仅包含高速公路的相应 OSM 文件。我想将OSM文件中的经纬度值转换为像素，并想在卫星图像上突出显示高速公路。

我尝试了 StackExchange 上介绍的几种方法。但是我得到每个经度和纬度值的负像素值和相同的像素值。有人可以解释一下，我错过了什么吗？

这是我使用 OTB 应用程序收集的图像信息。

这是我正在使用的代码。

from osgeo import gdal, osr
import numpy as np
import xml.etree.ElementTree as xml

src_filename = 'image.tif'
dst_filename = 'foo.tiff'


def readLongLat(path):
    lonlatList = []
    latlongtuple = ()
    root = xml.parse(path).getroot()
    for i in root:
        if i.tag == "node":
            latlong = []
            lat = float(i.attrib["lat"])
            long = float(i.attrib["lon"])
            latlong.append(lat)
            latlong.append(long)
            lonlatList.append(latlong)

    return lonlatList

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Get raster projection
epsg = 4269         # http://spatialreference.org/ref/sr-org/lambert_conformal_conic_2sp/
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)

# Make WGS84 lon lat coordinate system
world_sr = osr.SpatialReference()
world_sr.SetWellKnownGeogCS('WGS84')

transform = src_ds.GetGeoTransform()
gt = [transform[0],transform[1],0,transform[3],0,-transform[5]]

#Reading the osm file
lonlat = readLongLat("highways.osm")
# Transform lon lats into XY
coord_transform = osr.CoordinateTransformation(world_sr, srs)
newpoints = coord_transform.TransformPoints(lonlat) # list of XYZ tuples

# Make Inverse Geotransform  (try:except due to gdal version differences)
try:
    success, inverse_gt = gdal.InvGeoTransform(gt)
except:
    inverse_gt = gdal.InvGeoTransform(gt)

# [Note 1] Set pixel values
marker_array_r = np.array([[255]], dtype=np.uint8)
marker_array_g = np.array([[0]], dtype=np.uint8)
marker_array_b = np.array([[0]], dtype=np.uint8)
for x,y,z in newpoints:
    pix_x = int(inverse_gt[0] + inverse_gt[1] * x + inverse_gt[2] * y)
    pix_y = int(inverse_gt[3] + inverse_gt[4] * x + inverse_gt[5] * y)
    dst_ds.GetRasterBand(1).WriteArray(marker_array_r, pix_x, pix_y)
    dst_ds.GetRasterBand(2).WriteArray(marker_array_g, pix_x, pix_y)
    dst_ds.GetRasterBand(3).WriteArray(marker_array_b, pix_x, pix_y)

# Close files
dst_ds = None
src_ds = None

Answer 1

我最近尝试过使用xarray 模块。我认为xarray 是pandas 和numpy 之间的混合体，它允许您将信息存储为数组，但只需使用.sel 请求即可访问它。文档here.

更新：似乎需要安装rasterio 和xarray 才能使以下方法起作用。见link。

这是一种将 GeoTiff 文件转换为用户友好数组的更简单的方法。请参阅下面的示例：

import xarray as xr

ds = xr.open_rasterio("/path/to/image.tif")

# Insert your lat/lon/band below to extract corresponding pixel value
ds.sel(band=2, lat=19.9, lon=39.5, method='nearest').values
>>> [10.3]

这并不能直接回答您的问题，但可以帮助您确定我最近改用的不同（并且可能更简单）的方法。

注意：显然需要注意确保您的纬度/经度对与 GeoTiff 文件位于同一坐标系中，但我认为您无论如何都在处理它。

Answer 2

我可以使用 geoio 库来做到这一点。

import geoio  

img = geoio.GeoImage(src_filename)
pix_x, pix_y = img.proj_to_raster(lon,lat)