c/c++ 不使用GDAL读取geotiff全球高程数据

2025-04-13 23:33:09

背景

从tiff数据中按经纬度读取高程数据，要让依赖环境尽量简单，查阅资料GDAL依赖Boost，太重了,全球高程tif从www.ncei.noaa.gov/products/et… 下载

依赖库

libtiff libgeotiff libproj zlib libcurl sqlite

#include <geotiff.h>
#include <geotiffio.h>
#include <proj.h>
#include <tiffio.h>
#include <iostream>
#include "xtiffio.h"
#include "geo_normalize.h"
void MyTIFFErrorHandler(const char* module, const char* fmt, va_list args) {
    // 格式化错误消息
    char buffer[1024];
    vsnprintf(buffer, sizeof(buffer), fmt, args);
    // 输出到标准错误（或日志文件）
    std::cerr << "[TIFF ERROR] Module: " << (module ? module : "unknown")
              << ", Message: " << buffer << std::endl;
}
float ReadFloatTifGrid(const char *file,double lon, double lat ) {
    TIFFErrorHandler oldHandler = TIFFSetErrorHandler(MyTIFFErrorHandler);
    TIFF *tif = NULL;
    GTIF *gtif = NULL;
    tif = XTIFFOpen(file, "r");
    if (!tif) {
        return -99999 ;
    }
    gtif = GTIFNew(tif);
    if (!gtif) {
         GTIFFree(gtif);
        XTIFFClose(tif);
        return -99999  ;
    }

    unsigned short sampleFormat, samplesPerPixel, bitsPerSample;
    TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &samplesPerPixel);// 1: 无符号整数，2: 有符号整数，3: 浮点
    TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bitsPerSample);
    TIFFGetField(tif, TIFFTAG_SAMPLEFORMAT, &sampleFormat);

    if (sampleFormat != SAMPLEFORMAT_IEEEFP || bitsPerSample != 32 || samplesPerPixel != 1) {
        GTIFFree(gtif);
        XTIFFClose(tif);
        return -99999 ;
    }
    //图片大小
    int width, height;
    TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &width);
    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &height);

    short tiepointsize, pixscalesize;
    double* tiepoints;//[6];
    double* pixscale;//[3];
    if(!TIFFGetField(tif, TIFFTAG_GEOTIEPOINTS,  &tiepointsize,
                      &tiepoints)||
        !TIFFGetField(tif, TIFFTAG_GEOPIXELSCALE, &pixscalesize, &pixscale))
    {
        GTIFFree(gtif);
        TIFFClose(tif);
        return -99999 ;
    }
    // 3. 获取 GeoTIFF 定义
    GTIFDefn defn;
    if (!GTIFGetDefn(gtif, &defn)) {
        std::cerr << "Could not parse GeoTIFF definition." << std::endl;
        GTIFFree(gtif);
        TIFFClose(tif);
        return -99999 ;
    }
    // 4. 坐标转换（WGS84经纬度 -> 图片坐标系）
    char *target_crs = NULL;
    if (defn.Model == ModelTypeProjected && defn.PCS != KvUserDefined) {
        target_crs = (char *)malloc(16);
        snprintf(target_crs, 16, "EPSG:%d", defn.PCS);
    } else if (defn.Model == ModelTypeGeographic && defn.GCS != KvUserDefined) {
        target_crs = (char *)malloc(16);
        snprintf(target_crs, 16, "EPSG:%d", defn.GCS);
    } else {
        // 自定义坐标系处理（示例：UTM zone 50N）
        target_crs = "+proj=utm +zone=50 +north +datum=WGS84";
    }
    PJ_CONTEXT *ctx = proj_context_create();
    PJ *trans = proj_create_crs_to_crs(ctx, "EPSG:4326", target_crs, NULL);
    PJ_COORD proj = proj_trans(trans, PJ_FWD, proj_coord(lon, lat, 0, 0));

    // 5. 计算经纬度对应的像素坐标
    double A = tiepoints[3];  // 左上角地理坐标X
    double C = tiepoints[4];  // 左上角地理坐标Y
    double B = pixscale[0];   // 东西向像素分辨率
    double D = -pixscale[1];  // 南北向像素分辨率（取负号，因为Y轴向下）
    int x_pixel = (proj.xy.x - A) / B;
    int y_pixel = (proj.xy.y - C) / D;
     // 5. 读取数据
    double value=-99999;
    if (TIFFIsTiled(tif)) {//块存储
        // 处理分块图像
        uint32_t tile_width, tile_height;
        TIFFGetField(tif, TIFFTAG_TILEWIDTH, &tile_width);
        TIFFGetField(tif, TIFFTAG_TILELENGTH, &tile_height);

        unsigned char* buf = (unsigned char*)_TIFFmalloc(TIFFTileSize(tif));

       //unsigned char* buf1 = (unsigned char*)_TIFFmalloc(TIFFTileSize(tif));

        TIFFReadTile(tif, buf, x_pixel, y_pixel, 0, 0);
        // TIFFReadTile(tif, buf1,x_pixel+1, y_pixel, 0, 0);//即使偏移buf也是一样的 因为是按块读取
        if (bitsPerSample == 32) {
            float* floatData = (float*)buf;
            uint32_t offset_x = x_pixel % tile_width;
            uint32_t offset_y = y_pixel % tile_height;
            value = (floatData)[offset_y * tile_width + offset_x];

            // float* floatData1 = (float*)buf1;
            // float value1 = (floatData1)[offset_y * tile_width + offset_x];
            // int  a=0;
            // 处理32位浮点数据...
        } else if (bitsPerSample == 64) {
            double* doubleData = (double*)buf;
  
            // 处理64位浮点数据...
        }
        else if (bitsPerSample == 16) {
            uint16_t* shortData = (uint16_t*)buf;
        }
        //读取所有
        // for (uint32_t y = 0; y < height; y += tileHeight) {
        //     for (uint32_t x = 0; x < width; x += tileWidth) {
        //         TIFFReadTile(tif, buf, x, y, 0, 0);
        //         if(y>=3437&&x>18042)
        //         {
        //             int a=0;
        //         }
        //         if (bitsPerSample == 32) {
        //             float* floatData = (float*)buf;
        //             int a=0;
        //             // 处理32位浮点数据...
        //         } else if (bitsPerSample == 64) {
        //             double* doubleData = (double*)buf;
        //             int a=0;
        //             // 处理64位浮点数据...
        //         }
        //         else if (bitsPerSample == 16) {
        //             uint16_t* shortData = (uint16_t*)buf;

        //         }

        //     }
        // }
       _TIFFfree(buf);
    }
    else
    {
        int size=TIFFTileSize(tif);
        unsigned char* buf = (unsigned char*)_TIFFmalloc(TIFFScanlineSize(tif));
        for (long i = 0; i < height; i++) {
            if (TIFFReadScanline(tif, buf, (unsigned int)i, 0) == -1) {

                GTIFFree(gtif);
                XTIFFClose(tif);
                return value;
            }
        }
    }
    GTIFFree(gtif);
    XTIFFClose(tif);
    return value;
}

int main(int argc, char *argv[]) {

    float value=ReadFloatTifGrid("D:/Download/ETOPO_2022_v1_60s_N90W180_surface.tif",115.7,39.4);
    return 0;
}

注意

自己编译tiff库要开启zlib

总结

大模型不是那么好，生成代码乱七八糟，用于查查API还可以