如何为ASCII Grid格式地图数据配置正确的经纬度坐标轴？

我尝试用ASCII Grid和GML格式的地图数据绘图，希望x/y轴显示经纬度，但当前绘图轴是无法理解的默认值。代码如下：

# join ASCII tiles together to get area of interest
tile_list = glob.glob("data/01/input/OS/NY/NY[0-9][0-9].asc")
vrt = gdal.BuildVRT("merged.vrt", tile_list)
gdal.Translate("mergeDEMvrt.tif", vrt)
vrt = None

grid = Grid.from_raster('mergeDEMvrt.tif', data_name='dem')
dem = grid.read_raster('mergeDEMvrt.tif')

fig, ax = plt.subplots(figsize=(12,10))
fig.patch.set_alpha(0)

plt.imshow(dem, extent=grid.extent, cmap='terrain', zorder=1)
plt.colorbar(label='Elevation (m)')
plt.grid(zorder=0)
plt.title('Digital elevation map', size=14)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.tight_layout()

绘图轴格式不正确，需要显示常规经纬度数值。合并单个瓦片时是否需要叠加其他图层？我曾尝试另一个tif文件，其地理坐标能正常显示。也可获取shapefile、GML、GeoPackage或矢量瓦片格式的文件，同样遇到问题。

• 检查合并后TIFF的地理配准信息
  先确认合并后的mergeDEMvrt.tif是否带有正确的地理坐标系和投影信息，用GDAL命令行工具验证：

  gdalinfo mergeDEMvrt.tif
  

  若输出无Coordinate System相关内容，说明合并过程丢失了地理信息。需确保原始ASCII瓦片头部包含地理范围信息，若没有，需手动指定坐标系后再合并。

• 替换Grid类的地理范围读取逻辑
  部分Grid实现可能未正确解析TIFF的地理变换参数，改用GDAL直接读取：

  from osgeo import gdal
  
  ds = gdal.Open('mergeDEMvrt.tif')
  geotransform = ds.GetGeoTransform()
  # 计算实际地理范围
  x_min = geotransform[0]
  y_max = geotransform[3]
  x_max = x_min + geotransform[1] * ds.RasterXSize
  y_min = y_max + geotransform[5] * ds.RasterYSize
  extent = (x_min, x_max, y_min, y_max)
  ds = None
  
  # 绘图时替换原extent
  plt.imshow(dem, extent=extent, cmap='terrain', zorder=1)
  

• 投影坐标系转经纬度
  若原始数据是投影坐标系（如UTM），需转换为WGS84经纬度，用pyproj实现：

  import pyproj
  from matplotlib.ticker import FormatStrFormatter
  
  # 替换为原始数据的EPSG编码（从gdalinfo获取）
  src_proj = pyproj.CRS.from_epsg(32618)
  dst_proj = pyproj.CRS.from_epsg(4326)  # WGS84经纬度
  transformer = pyproj.Transformer.from_crs(src_proj, dst_proj, always_xy=True)
  
  # 转换范围角点
  x_min_lon, y_min_lat = transformer.transform(x_min, y_min)
  x_max_lon, y_max_lat = transformer.transform(x_max, y_max)
  lon_lat_extent = (x_min_lon, x_max_lon, y_min_lat, y_max_lat)
  
  # 绘图并格式化坐标轴
  plt.imshow(dem, extent=lon_lat_extent, cmap='terrain', zorder=1)
  ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f°'))
  ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f°'))
  

• 合并时保留地理信息
  合并瓦片无需叠加额外图层，关键是确保合并过程保留地理配准。构建VRT时可指定投影：

  vrt_options = gdal.BuildVRTOptions(outputSRS='EPSG:4326')  # 替换为实际坐标系
  vrt = gdal.BuildVRT("merged.vrt", tile_list, options=vrt_options)
  

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