当前位置：首页 > news >正文

从稀疏数据（CSV）创建非常大的 GeoTIFF（和 WMS）

news 2025/7/27 22:55:33

总体流程

首先，数据是基于点的：一个经纬度点的数据值。我们生成这些数据并输出为 CSV 文件。根据分辨率和我们评分的区域，即使数据非常稀疏，CSV 文件中也会有数亿条记录。然后，您需要经过多个步骤，才能将 CSV 文件转换为 GeoTIFF 格式，再转换为 WMS 地图图块：

生成数据，生成一个包含纬度、经度和值（也可以包含其他列）的 CSV 文件。
将 CSV 转换为XYZ 文件，该文件对列进行排序，并且每个点之间的间距必须相等（而 CSV 中的经纬度点可以按任意顺序排列）。“z”值是您的数据值，通常是海拔、温度等。在我们的例子中，它通常是NatureScore®。
从 XYZ 创建一个 GeoTIFF 文件，并确保其中存在一个无数据值（由于数据稀疏，会有很多像素没有数据值）。GeoTIFF 只是一个 TIFF 文件，但其中每个像素都有一个 GIS 坐标值（在我们的例子中是经纬度 (EPSG 4326)，但也可能是基于 CRS（坐标参考系）的其他数字）。
根据您的需要为 GeoTIFF 着色（即将您的数据值转换为颜色，因为 WMS 是基于图像的服务，并且您无法指定如何将值映射到颜色）。
从 GeoTIFF 生成 WMS 地图图块图像。
使用地图图块创建您的 WMS 服务器。

问题

这个过程的大部分内容都很常见，并且有详尽的文档。然而，我在第 3 步（创建 GeoTIFF）遇到了一个问题，因为 XYZ 数据非常庞大且稀疏。下面是可以使用的命令：

gdal_translate -a_srs EPSG:4326 -of GTiff -ot Byte ${xyzfile} ${geotiff}

运行不会出错，但生成的 GeoTIFF 文件为空或不正确。我找不到之前遇到的那个关于处理gdal_translate非常大的稀疏 CSV 文件的问题，但这似乎是一个已知问题。

为了解决这个问题，我创建了一个相当简单的 Python 程序，用零值填充所有缺失的坐标，这解决了问题，但也创建了一个非常庞大的 XYZ 文件（例如，仅对 3 个城市的数据进行小规模测试，就生成了一个 45GB 的文件，因此你可以想象整个美国的文件会有多大）。不过，我反复进行了迭代，借助 Claude AI 的帮助，将原始解决方案调整为更简单的解决方案，并且不会在大数据量下崩溃。原始方案使用了 scipy 和 rasterio，它会在完整数据集上崩溃。新的解决方案仅使用了 numpy、pandas 和 gdal。我还跳过了输出更新后的 XYZ 文件的步骤。但真正的关键在于：

numpy.full((rows, cols), NODATA)

我的 NODATA 值为 255。其他值均为 1-10。然后程序会写入 GeoTIFF 文件。写入这些 GeoTIFF 文件时，务必使用 LZW 压缩，因为尤其是在处理稀疏数据时，它会显著减小文件大小。在我的例子中，未压缩时 GeoTIFF 文件大小为 74GB，而压缩后则下降到 187MB！

一旦我有了这个 GeoTIFF，我就可以顺利完成其余的过程……

全部细节/流程

我会按步骤分解这个过程。当然，你可以自动化整个过程（例如，我有一个shell脚本，它接收CSV文件并生成中间文件（XYZ、GeoTIFF），最终生成地图图块），但你的需求可能会有所不同，等等。

首先，您需要准备好 CSV 或其他 GIS 源数据集。我将使用 CSV 进行演示，但实际上，无论您拥有什么数据（假设是栅格数据），您只需要将其转换为 XYZ 格式即可。您可以通过 QGIS 将数据集导出为 XYZ 格式等方式来实现。也就是说，根据您的数据，此步骤可能需要也可能不需要。

CSV 到 XYZ

以下是将 CSV 文件转换为 XYZ 格式的 Shell 脚本的基础知识。正如您将看到的，我的源数据包含纬度、经度和其他几列，但第 8 列是需要保留的数据 (z)。该脚本会将数据裁剪为仅保留这几列，对 y、x 列进行排序，并将列标题重命名为 y,x,z（纬度为 y，经度为 x，数据为 z）。注意，即使是“XYZ”文件，您的列顺序也可能不同：

inputcsv=$1
xyzfile=${inputcsv//.csv/.xyz}# Sort our data, and trim it to be just XYZ columns
# latitude: column 1
# longitude: column 2
# our data: column 8
cat ${inputcsv} | sort -n -t ',' -k1 -k2 | cut -d, -f1-2,8 > ${xyzfile}# Rename columns to be YXZ:
sed -i "" -e 's/latitude/y/' ${xyzfile}
sed -i "" -e 's/longitude/x/' ${xyzfile}
sed -i "" -e 's/Data/z/' ${xyzfile}

XYZ 到 GeoTIFF

接下来，将 XYZ 数据转换为单波段 GeoTIFF。我不会在这里列出所有内容，但这里有一个 Python 程序的要点。你可能需要根据自己的需求调整常量NODATA。此外，如果你不太熟悉 Python，可以通过安装 Python 3 来设置环境，然后创建一个虚拟环境来指定版本和依赖项，例如：

python3 -m venv .
source bin/activatepip install pandas
pip install numpy
pip install gdal
pip install osr

注意，你需要安装 GDAL（3.11 或更高版本）。在 Mac 上，你可以使用来完成此操作brew install gdal。

然后执行以下操作：

python xyz_to_geotiff.py your-data.xyz your-data.tiff

请注意，该脚本还提供了指定分辨率和 SRID 的选项。它会尝试为您确定分辨率——我知道自己所需的分辨率，因此指定了分辨率以跳过该步骤。该脚本确实会写入压缩的 GeoTIFF 文件，但我最初并没有这样做，因此统计数据显示初始文件大小为 74GB。

为 GeoTIFF 着色

为了给 GeoTIFF 文件着色，使其更易于查看，并可能具有更清晰的含义，您可以使用gdaldem将数据值映射到颜色的颜色文件。颜色文件是一种简单的映射，其格式为数据值、空格和颜色值。例如，以下是 NatureQuant 红绿调色板和 255/无数据值（黑色）：

10  204:51:51
9  221:68:68
8  238:102:51
7  255:153:0
6  255:204:0
5  255:221:0
4  221:221:0
3  153:204:0
2  0:153:51
1  0:102:51
255 0:0:0

如你所见，顺序无关紧要，它只是提供了一张用于着色的地图。因此，要生成 GeoTIFF 的彩色版本，请执行以下操作：

gdaldem color-relief -co COMPRESS=LZW \your-data.tiff \your-colors.txt  \ your-data-colored.tiff

这将输出一个新的 GeoTIFF，其中包含 3 个波段（RGB）的数据，每个像素现在都具有与其值对应的颜色your-colors.txt。以下是结果示例：

没有任何

这就是创建 WMS 所需的全部内容，但让我们再介绍一个可选步骤……

在 GeoTIFF 中使用透明度（针对无数据区域）

上面的图片没问题，但是，根据你的需求，也许一个版本能使所有无数据区域都透明更有意义？如果是这样，你可以按照以下步骤操作：

gdal_translate -a_nodata 0 your-data-colored.tiff your-data-nodata0.tiffgdal_translate -of GTiff -co "ALPHA=YES" your-data-nodata0.tiff your-data-with-alpha.tiff

这会添加第四个波段（一个 alpha 蒙版）并将您的无数据值映射到它，从而生成一个新的 GeoTIFF，其中所有黑色区域现在都是 100% 透明的。

为 WMS 创建地图图块

最后一步是为 WMS 服务器创建地图图块图像。这很简单：

mkdir map_tilesgdal2tiles --xyz -z 4-13 --srcnodata=0 --resume your-data-colored.tiff map_tiles

在本例中，我们指定创建 4-13 级缩放级别，然后指定--resume参数。我发现这需要运行几次才能完成，gdal2tiles有时会崩溃。幸好有了 resume 参数，它基本上会从上次中断的地方继续执行（并且会跳过生成已存在的文件——它不会对现有文件进行健全性检查，而是纯粹根据文件的存在情况来决定是否应该处理该文件）。

对于我的数据，这在 M3 Macbook Pro 上花费了大约 10 个小时，并产生了近460 万个文件（18GB）。

你可以使用程序在 map_tiles 文件夹中生成的文件来验证效果leaflet.html（只需在浏览器中打开该 HTML 文件即可），你应该可以看到完整的地图。以下是我们部分数据的示例：

没有任何