将Numpy数组转为带地理参考框的GeoTIFF：解决错位转置问题

你的问题出在地理变换参数顺序错误和数组存储方向与GeoTIFF规范不匹配两个核心点上，以下是具体修正步骤：

问题根源 #

1. 位置偏移：rasterio.transform.from_bounds的参数顺序是west, south, east, north，但你传入的extent是[xedges[0], xedges[-1], yedges[0], yedges[-1]]（即west, east, south, north），参数顺序完全错误，导致地理范围偏移。
2. 图像转置/颠倒：matplotlib.imshow用origin='lower'时，数组是从南到北（下到上）排列，但GeoTIFF要求行顺序是从北到南（上到下），两者方向相反，导致图像翻转。

修正后的完整代码 #

1. 调整保存GeoTIFF的核心代码 #

import numpy as np
import rasterio
from rasterio.transform import from_bounds

# 假设img和extent来自myplot函数的返回（以某个sigma值为例）
img, extent = myplot(x, y, weights, s=32)

# 修正1：翻转数组y轴，匹配GeoTIFF北->南的行顺序
img_geotiff = np.flipud(img)

# 修正2：使用正确的from_bounds参数顺序：west, south, east, north
# 同时用图像实际维度代替硬编码的10000，避免维度不匹配
transform = from_bounds(
    extent[0],    # west
    extent[2],    # south
    extent[1],    # east
    extent[3],    # north
    img_geotiff.shape[1],  # 图像宽度（x方向像素数）
    img_geotiff.shape[0]   # 图像高度（y方向像素数）
)

# 修正3：保存GeoTIFF，用nan作为无数据值（避免与有效0密度混淆）
with rio.open(f'data_sigma_{s}.tif', 'w', 
              driver='GTiff', 
              height=img_geotiff.shape[0], 
              width=img_geotiff.shape[1], 
              count=1, 
              dtype='float64', 
              nodata=np.nan,
              crs='EPSG:32632', 
              transform=transform) as dst:
    dst.write(img_geotiff, 1)

2. 循环内批量保存的修改 #

如果要在原来的循环中批量生成GeoTIFF，可以把保存逻辑嵌入循环：

for ax, s in zip(axs.flatten(), sigmas):
    if s == 0:
        ax.plot(x, y, weights, 'k.', markersize=5)
        ax.set_title("Scatter plot")
        plt.savefig('export_'+str(s)+'.png', dpi=150, bbox_inches='tight')
    else:
        img, extent = myplot(x, y, weights, s)
        ax.imshow(img, extent=extent, origin='lower', cmap=cm.jet)
        ax.set_title(f"Smoothing with $\sigma$ = {s}")
        plt.savefig(f'export_{s}.png', dpi=150, bbox_inches='tight')
        
        # 新增：保存为GeoTIFF
        img_geotiff = np.flipud(img)
        transform = from_bounds(extent[0], extent[2], extent[1], extent[3], img.shape[1], img.shape[0])
        with rio.open(f'data_sigma_{s}.tif', 'w', 
                      driver='GTiff', 
                      height=img.shape[0], 
                      width=img.shape[1], 
                      count=1, 
                      dtype='float64', 
                      nodata=np.nan,
                      crs='EPSG:32632', 
                      transform=transform) as dst:
            dst.write(img_geotiff, 1)

关键修改说明 #

• np.flipud(img)：将数组沿y轴上下翻转，把matplotlib中南->北的顺序转为GeoTIFF要求的北->南顺序，解决图像颠倒问题。
• from_bounds参数修正：把原代码中的extent[1]（east）和extent[2]（south）交换位置，确保地理范围正确映射。
• 用img.shape代替硬编码的10000：如果后续修改bins参数，图像维度会变化，硬编码值会导致地理变换分辨率不匹配。
• 改用np.nan作为nodata：0可能是有效的住宅密度值，用nan更准确区分无数据区域。

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