转换问题：大型 TIFF 文件到 XYZ 格式

Question

我正在使用名为 CHELSAcruts_tmax_4_1981_V.1.0.tif（大小：97M）的 TIFF 文件。 我想先将此文件转换为 XYZ 格式，然后再转换为 CSV。

尝试与挑战 在线转换器限制：最初，我尝试使用在线转换器，但文件大小超出了允许的限制。

用于 XYZ 转换的 Python 脚本： 我找到了以下用于转换为 XYZ 格式的 Python 代码片段：

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

但是，这种方法有两个主要问题：

运行时间过长。 它生成的文件大于我的虚拟机的可用磁盘空间（超过 100 GB）。 零值的观察：经过仔细检查，我注意到该文件包含许多零值，例如：

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

文档：https://gdal.org/api/python/osgeo.gdal.html 尝试过滤零值：为了解决这个问题，我尝试使用以下代码过滤这些值，但不成功：

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()

我尝试使用以下代码清除零值：

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)

但它也永远有效并且没有任何结果

https://drive.google.com/file/d/1DPpxNFgcE3d4W-7AemBN_zTtylibI9Na/view?usp=sharing

鉴于尺寸较大，如何有效地将其转换为 XYZ，然后转换为 CSV 格式，同时管理尺寸和零值问题？

Answer 1

出于多种原因，这样做对我来说没有任何意义。 XYZ/CSV 格式将存储每个有效像素的坐标（作为字符串！），这对于像这样的常规网格来说是完全多余的。可以使用例如地理变换（6 个值）和文件尺寸（2 个值）来检索完全相同的信息。

将其存储为字符串，未压缩会完全耗尽数据量。使用不同的格式（例如（地理）镶木地板）已经可以缓解很多问题。

无论如何，下面的代码片段展示了如何完成它。这会一次性完成所有操作，因此您需要内存来执行此操作。当然可以调整代码并以较小的块读取数据，然后将它们增量写入输出文件。

在这种情况下使用 Pandas，为了方便起见，可能还会在此基础上引入它自己的开销。对于像这样的简单文件格式，手动执行也是可行的。

from osgeo import gdal
import numpy as np
import pandas as pd

读取栅格数据和元数据

fn = "CHELSAcruts_tmax_4_1981_V.1.0.tif"

ds = gdal.Open(fn)

gt = ds.GetGeoTransform()

data = ds.ReadAsArray()
ys, xs = data.shape

ds = None

过滤节点数据并创建坐标

valid = data != -32768
valid.sum(), valid.size

ulx, uly = gdal.ApplyGeoTransform(gt, 0.5, 0.5)
lrx, lry = gdal.ApplyGeoTransform(gt, xs-0.5, ys-0.5)

lats, lons = np.mgrid[uly:lry:gt[5], ulx:lrx+gt[1]:gt[1]]

# reduce memory a little...
lons = lons.astype(np.float32)
lats = lats.astype(np.float32)

写入 csv

df = pd.DataFrame(dict(
    X=lons[valid],
    Y=lats[valid],
    Z=data[valid],
))

df.to_csv("CHELSAcruts_tmax_4_1981_V.1.0.csv", index=False)

如前所述，我会考虑另一种方法，至少将数据存储为二进制、压缩的且没有冗余坐标。