大体积TIFF转XYZ及CSV格式的高效处理方案求助

我正在处理CHELSAcruts_tmax_4_1981_V.1.0.tif（大小97M），需要先转XYZ格式再转CSV格式，但遇到了以下问题：

尝试过的方法与问题 #

• 在线转换器：文件大小超出平台上限，无法使用。

• GDAL Python脚本转XYZ：
  用以下代码转换时，存在两个严重问题：

  from osgeo import gdal
  ds = gdal.Open("CHELSAcruts_tmax_4_1981_V.1.0.tif")
  ds.GetGeoTransform()
  xyz = gdal.Translate("dem.xyz", ds)
  

  • 运行时间极长；
  • 生成的XYZ文件超过VM可用磁盘空间（超100GB）。

• 无数据值问题：文件包含大量-32768无数据值，示例：

  X Y Z
  -179.99597222220001 83.9956937485000168 -32768
  -179.987638888900022 83.9956937485000168 -32768
  

  尝试通过GDAL过滤无数据值但失败：

  band = ds.GetRasterBand(1)
  band.SetNoDataValue(-32768)
  gdal.Translate("dem.xyz", ds, noData=-32768, creationOptions=["ADD_HEADER_LINE=YES"])
  

• 查看数据规模：通过代码确认图像数组形状为(20880, 43200)，这是文件膨胀的核心原因：

  from osgeo import gdal
  ds = gdal.Open('mytif.tif', gdal.GA_ReadOnly)
  rb = ds.GetRasterBand(1)
  img_array = rb.ReadAsArray()
  

• 尝试替换无数据值：用gdal_calc替换-32768为NaN，但代码运行无响应，未生成结果：

  import gdal_calc
  import numpy as np
  
  original_tiff = "/content/drive/My Drive/Kalmia/CHELSAcruts_tmax_4_1981_V.1.0.tif"
  modified_tiff = "/content/drive/My Drive/Kalmia/modified_tiff.tif"
  
  # Use gdal_calc.py to replace -32768 with NaN
  gdal_calc.Calc("A*(A!=-32768) + nan*(A==-32768)", A=original_tiff, outfile=modified_tiff, NoDataValue=np.nan)
  

现在需要高效完成TIFF→XYZ→CSV的转换，同时解决文件过大和无数据值过滤的问题。

1. 直接过滤无数据值生成XYZ（避免中间文件膨胀） #

问题出在gdal.Translate的参数用法错误，正确做法是明确指定源无数据值并让GDAL跳过这些像素，而非写入后再处理。

使用优化后的GDAL Python代码：

from osgeo import gdal

# 打开TIFF文件
ds = gdal.Open("CHELSAcruts_tmax_4_1981_V.1.0.tif", gdal.GA_ReadOnly)

# 设置无数据值，确保GDAL识别
band = ds.GetRasterBand(1)
band.SetNoDataValue(-32768)

# 转换时跳过无数据像素，生成过滤后的XYZ
xyz_options = gdal.TranslateOptions(
    format="XYZ",
    srcNodata=-32768,  # 标记源数据的无数据值
    dstNodata=None,    # 目标不保留无数据值，直接跳过
    creationOptions=["ADD_HEADER_LINE=YES"]
)

gdal.Translate("filtered_dem.xyz", ds, options=xyz_options)

核心作用：srcNodata=-32768让GDAL直接忽略无数据像素，不会写入到XYZ文件中，大幅缩减文件体积，同时减少不必要的计算。

2. 直接转CSV（跳过XYZ中间步骤，更高效） #

如果最终目标是CSV，完全可以跳过XYZ格式，直接从TIFF生成CSV，进一步节省时间和空间：

import numpy as np
from osgeo import gdal

ds = gdal.Open("CHELSAcruts_tmax_4_1981_V.1.0.tif", gdal.GA_ReadOnly)
geo_transform = ds.GetGeoTransform()
band = ds.GetRasterBand(1)
band.SetNoDataValue(-32768)

# 读取数据数组
data = band.ReadAsArray()

# 获取地理坐标参数
x_min = geo_transform[0]
x_res = geo_transform[1]
y_max = geo_transform[3]
y_res = geo_transform[5]

# 生成所有像素的X/Y坐标网格
x_coords = np.arange(x_min, x_min + data.shape[1] * x_res, x_res)
y_coords = np.arange(y_max, y_max + data.shape[0] * y_res, y_res)

# 过滤无数据值，提取有效数据和对应坐标
mask = data != -32768
flat_data = data[mask]
flat_x = x_coords[np.tile(mask, (1,1))].reshape(data.shape)[mask]
flat_y = y_coords.repeat(data.shape[1])[mask.ravel()]

# 组合成CSV格式并写入
csv_data = np.column_stack((flat_x, flat_y, flat_data))
np.savetxt("output.csv", csv_data, delimiter=",", header="X,Y,Z", comments='')

优势：

• 无需生成庞大的XYZ中间文件；
• 利用numpy向量化操作，比逐像素处理快数倍；
• 直接过滤无数据值，生成的CSV仅包含有效数据。

3. 命令行GDAL工具（更快，适合超大文件） #

如果Python脚本仍有性能瓶颈，直接使用GDAL命令行工具，效率更高：

# 生成过滤后的XYZ文件
gdal_translate -of XYZ -srcnodata -32768 -dstnodata "" -co ADD_HEADER_LINE=YES CHELSAcruts_tmax_4_1981_V.1.0.tif filtered_dem.xyz

# 将XYZ转CSV（替换空格为逗号）
sed 's/ /,/g' filtered_dem.xyz > output.csv

关键参数：-dstnodata ""确保无数据像素不被写入，直接跳过，是控制文件体积的核心。

内存优化提示 #

如果VM内存不足，避免一次性加载整个数据集，可分块读取处理：

block_size = band.GetBlockSize()
# 逐块读取并处理
for y in range(0, data.shape[0], block_size[1]):
    for x in range(0, data.shape[1], block_size[0]):
        block = band.ReadAsArray(x, y, block_size[0], block_size[1])
        # 处理当前块的过滤和坐标生成，追加写入文件

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