坐标转像素栅格越界求助：从SUMO .net.xml提取坐标匹配DGM1.tif高程

问题背景 #

从SUMO的.net.xml文件提取Shape属性中的米单位x/y坐标，需要从DGM1文件夹的.tif文件中获取对应高程。目前坐标提取成功，但转换为像素坐标时触发"raster is out of bounds"错误。

栅格数据参数 #

• 像素参数：宽度1.0，旋转参数1/2均为0.0，高度-1.0
• 左上角地理坐标：东坐标462000.5，北坐标5548999.5
• 数据集属性：宽1000、高1000、单波段、CRS为EPSG:25832
• 边界框：left=462000.0，bottom=5548000.0，right=463000.0，top=5549000.0
• 数据类型：float32

当前错误的转换计算 #

用户采用的转换公式及示例计算：
X坐标转换：

X offset = X coordinate − 左上角东坐标

Y坐标转换：

Y offset = 左上角北坐标 − Y coordinate

以X=312.37、Y=1176.31为例：

X offset = 312.37−462000.5= −461688.13
Y offset =5548999.5−1176.31= 5547823.19

像素X坐标 = X offset / 像素宽度 = −461688.13/1.0 = −461688.13
像素Y坐标 = Y offset / |像素高度| = 5547823.19/1.0 = 5547823.19

计算结果远超出栅格0-999的有效范围，说明转换逻辑存在根本性错误。

现有Python代码 #

# Read the TIF file for elevation data
with rasterio.open(tif_file) as src:
    # Read elevation data
    elev_data = src.read(1)  # Assuming single-band elevation data

    # Read the TFW file for georeferencing information
    with open(tfw_file, 'r') as tfw:
        lines = tfw.readlines()
        # Extract georeferencing parameters from the TFW file
        pixel_width = float(lines[0])
        pixel_height = float(lines[3])
        x_coord = float(lines[4])
        y_coord = float(lines[5])

        # Compute the x and y coordinates in the raster
        col = int((x - x_coord) / pixel_width)
        row = int((y_coord - y) / pixel_height)
       # print (col, row)

        # Check if the coordinates are within the raster bounds
        if 0 <= row < src.height and 0 <= col < src.width:
            elevation = elev_data[row, col]
            return elevation

return None  # Return None if coordinates are not found in any file

问题根源与解决方案 #

问题根源 #

1. 坐标坐标系不匹配：SUMO提取的x/y大概率是模拟区域相对坐标（SUMO默认笛卡尔坐标系，原点在模拟区域内部），而栅格使用的是EPSG:25832大地坐标（UTM 32N的东/北向米单位），直接用相对坐标转换必然越界。
2. 手动转换公式误差：即使坐标单位匹配，手动计算未考虑栅格实际边界的起始偏移，且未利用rasterio内置的地理变换处理能力，容易出错。

修正步骤 #

步骤1：统一坐标系统

检查SUMO的.net.xml是否包含<location>标签，确认坐标系信息：

<location netOffset="0.00,0.00" convBoundary="0.00,0.00,1000.00,1000.00" origBoundary="-10000000000.00,-10000000000.00,10000000000.00,10000000000.00" projParameter="+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"/>

• 若存在projParameter且与栅格CRS一致，Shape中的x/y即为大地坐标；
• 若为相对坐标，需将提取的x/y加上netOffset的对应值，转换为大地坐标。

步骤2：使用rasterio内置转换（推荐）

无需手动读取TFW文件，直接用rasterio处理地理变换，避免计算错误：

with rasterio.open(tif_file) as src:
    # 先判断坐标是否在栅格边界内
    bbox = src.bounds
    if not (bbox.left <= x <= bbox.right and bbox.bottom <= y <= bbox.top):
        return None
    # 自动转换大地坐标到像素坐标
    try:
        row, col = src.index(x, y)
    except ValueError:
        return None
    # 读取对应高程
    return src.read(1)[row, col]

步骤3：手动修正转换公式（若需自定义）

当确认坐标为EPSG:25832大地坐标时，修正计算公式：

# 基于栅格实际边界计算，而非左上角偏移
col = int((x - 462000.0) / pixel_width)
# 像素高度为-1.0，直接用(上边界-Y)即可得到行号
row = int((5549000.0 - y) / abs(pixel_height))

内容的提问来源于stack exchange，提问作者Sumanth Rc