创建GeoTIFF的DataFrame值与提取值不匹配的原因及解决方法

问题现象 #

使用NumPy、Pandas、Rasterio创建GeoTIFF后，提取的像素值与原始DataFrame中的值完全不符，部分提取值为nodata（0.0），其余值也与原数据对应不上。

原因分析 #

1. 栅格行顺序颠倒 #

GeoTIFF规范要求栅格数据按**从上到下（纬度从高到低）**的顺序存储，但你的代码中：

• 创建的lat_mesh是从低纬度（39.0）到高纬度（40.9）排列
• 转换为z_array后直接写入GeoTIFF，导致低纬度数据被放到了GeoTIFF的顶部行，高纬度数据放到底部行，行顺序完全反转。

2. 采样点位置错误 #

你创建的采样点是经纬度网格的边界节点（如39.0, -76.0），但GeoTIFF的每个像素是一个区间（例如最下方像素的纬度范围是39.0~39.1），边界节点不在任何像素的有效范围内，因此Rasterio返回nodata值。而原始DataFrame中的z值对应的是像素中心的坐标，而非边界节点。

解决办法 #

修改GeoTIFF创建代码 #

在写入前将z_array沿行方向上下翻转，匹配GeoTIFF的存储顺序：

import numpy as np
import pandas as pd
import rasterio as rio

# 定义经纬度范围与分辨率
lat_range = (39, 41)
lon_range = (-76, -73)
resolution = 0.1

# 创建经纬度数组
lats = np.arange(lat_range[0], lat_range[1], resolution)
lons = np.arange(lon_range[0], lon_range[1], resolution)

# 生成网格并扁平化
lon_mesh, lat_mesh = np.meshgrid(lons, lats)
lat_values = lat_mesh.flatten()
lon_values = lon_mesh.flatten()

# 生成随机z值
z_values = np.random.rand(len(lat_values))
data = pd.DataFrame({'lat': lat_values, 'lon': lon_values, 'z': z_values})

# 定义GeoTIFF参数
width = len(lons)
height = len(lats)
transform = rio.transform.from_bounds(lon_range[0], lat_range[0], lon_range[1], lat_range[1], width, height)
crs = rio.crs.CRS.from_epsg(4326)

# 写入GeoTIFF：添加上下翻转步骤
with rio.open("output.tif", "w", driver="GTiff", width=width, height=height, count=1, dtype=np.float32, nodata=0, transform=transform, crs=crs) as dst:
    z_array = z_values.reshape((height, width))
    z_array = np.flipud(z_array)  # 关键：上下翻转数组，匹配GeoTIFF存储顺序
    dst.write(z_array, 1)

修改提取代码 #

使用像素中心坐标进行采样，并处理采样结果的格式：

import rasterio as rio
import pandas as pd
import geopandas as gpd
import numpy as np

with rio.open(r"E:\Machine_Learning_for_Himalaya_IEEE_GRSL\plots\output.tif") as src:
    print("meta:", src.meta)
    lat_range = (39, 41)
    lon_range = (-76, -73)
    resolution = 0.1

    # 生成像素中心的经纬度数组（而非边界节点）
    lats = np.arange(lat_range[0] + resolution/2, lat_range[1], resolution)
    lons = np.arange(lon_range[0] + resolution/2, lon_range[1], resolution)
    lon_mesh, lat_mesh = np.meshgrid(lons, lats)
    
    # 创建采样点并提取值
    geometry = gpd.points_from_xy(lon_mesh.ravel(), lat_mesh.ravel())
    points = gpd.GeoDataFrame(geometry=geometry)
    # 提取值并转换为单个数值（而非数组）
    zs = [val[0] for val in src.sample(zip(points.geometry.x, points.geometry.y))]

    # 生成结果DataFrame
    df = pd.DataFrame({'Latitude': lat_mesh.ravel(),
                       'Longitude': lon_mesh.ravel(),
                       'Value': zs})

print(df)

验证结果 #

修改后，提取的Value列将与原始DataFrame中的z列一一对应（需匹配像素中心坐标的行顺序），不会再出现nodata值，且数值完全一致。

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