使用 LibTIFF 读取和计算 GeoTIFF 高程数据

本博客文章提供了一个 C 代码示例，演示了如何使用 LibTIFF 库读取 GeoTIFF 文件并计算高程统计数据。

代码

#include <stdio.h>
#include <stdlib.h>
#include 'tiffio.h'

#define ROWS 327
#define COLS 486
#define PIXEL_SIZE 5
#define TILE_SIZE 64

int main(int argc, char* argv[]) {
    TIFF* tif;
    float* data;
    uint32_t width, height;
    uint16_t bits_per_sample, samples_per_pixel;
    tsize_t tile_size;
    float min_elevation, max_elevation;
    unsigned long long total_fill_volume = 0, total_cut_volume = 0;
    int i, j, tile_row, tile_col;

    if (argc < 2) {
        fprintf(stderr, 'Usage: %s <input file>
', argv[0]);
        exit(1);
    }

    tif = TIFFOpen(argv[1], 'r');
    if (!tif) {
        fprintf(stderr, 'Error opening TIFF file
');
        exit(1);
    }

    TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &width);
    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &height);
    TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bits_per_sample);
    TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &samples_per_pixel);

    if (width != COLS || height != ROWS || bits_per_sample != 32 || samples_per_pixel != 1) {
        fprintf(stderr, 'Invalid TIFF format
');
        exit(1);
    }

    tile_size = TIFFTileSize(tif);
    data = (float*)_TIFFmalloc(tile_size);
    if (!data) {
        fprintf(stderr, 'Error allocating memory
');
        exit(1);
    }

    if (!TIFFIsTiled(tif)) {
        fprintf(stderr, 'Can not read tiles from a non-tiled image.
');
        exit(1);
    }

    // 初始化最小高程为第一个读取的高程值
    TIFFReadTile(tif, data, 0, 0, 0, 0);
    min_elevation = data[0]; 
    max_elevation = data[0];

    for (tile_row = 0; tile_row < height; tile_row += TILE_SIZE) {
        for (tile_col = 0; tile_col < width; tile_col += TILE_SIZE) {
            TIFFReadTile(tif, data, tile_col, tile_row, 0, 0);
            for (i = 0; i < TILE_SIZE && tile_row + i < height; i++) {
                for (j = 0; j < TILE_SIZE && tile_col + j < width; j++) {
                    float elevation = data[i * TILE_SIZE + j];
                    if (elevation > max_elevation) {
                        max_elevation = elevation;
                    }
                    if (elevation < min_elevation) {
                        min_elevation = elevation;
                    }
                    if (elevation > 0) {
                        total_fill_volume += (PIXEL_SIZE * PIXEL_SIZE * elevation);
                    }
                    else {
                        total_cut_volume += (PIXEL_SIZE * PIXEL_SIZE * (-elevation));
                    }
                }
            }
        }
    }

    printf('Min elevation: %.3f meters\n', min_elevation);
    printf('Max elevation: %.3f meters\n', max_elevation);
    printf('Total fill volume: %.3f cubic kilometers\n', total_fill_volume / 1e9);
    printf('Total cut volume: %.3f cubic kilometers\n', total_cut_volume / 1e9);

    _TIFFfree(data);
    TIFFClose(tif);

    return 0;
}

解释

这段代码首先打开一个 GeoTIFF 文件并读取其元数据。然后，它遍历图像中的每个像素，计算最小和最大高程值，以及填充体积和切割体积。

最小高程初始化：

为了解决最小高程计算错误的问题，代码现在首先读取第一个图块的数据，并将 min_elevation 初始化为第一个读取的高程值。这样做可以确保 min_elevation 正确初始化，即使数据中存在负高程值。

其他计算:

填充体积是指高程高于 0 的区域的体积。
切割体积是指高程低于 0 的区域的体积。

结论

这段代码提供了一种简单的方法来使用 LibTIFF 库读取 GeoTIFF 文件并计算高程统计数据。它还展示了如何通过正确初始化最小高程来解决可能出现的错误。