c++ - 使用 GDAL 在 FAA Sectional 上查找像素位置

Question

我正在尝试将纬度和经度位置转换为 GEOTIFF FAA 剖面图中的像素位置 (x,y)。作为引用，这里是我正在使用的西雅图部分文件的一些 GDALInfo 输出。

    Driver: GTiff/GeoTIFF
    Files: d:\programming\gps\Seattle SEC 94.tif
    Size is 17869, 12333
    Coordinate System is:
    PROJCS["Lambert Conformal Conic",
        GEOGCS["NAD83",
            DATUM["North_American_Datum_1983",
                SPHEROID["GRS 1980",6378137,298.2572221010042,
                   AUTHORITY["EPSG","7019"]],
                AUTHORITY["EPSG","6269"]],
            PRIMEM["Greenwich",0],
            UNIT["degree",0.0174532925199433],
            AUTHORITY["EPSG","4269"]],
        PROJECTION["Lambert_Conformal_Conic_2SP"],
        PARAMETER["standard_parallel_1",46.66666666666666],
        PARAMETER["standard_parallel_2",41.33333333333334],
        PARAMETER["latitude_of_origin",46.75],
        PARAMETER["central_meridian",-121],
        PARAMETER["false_easting",0],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]]]
    Origin = (-409375.329186913440000,261032.079643933710000)
    Pixel Size = (42.334884781508684,-42.335597515618694)
    Image Structure Metadata:
      COMPRESSION=LZW
     INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  ( -409375.329,  261032.080) (126d34'48.90"W, 48d58'12.69"N)
    Lower Left  ( -409375.329, -261092.845) (126d 8'14.67"W, 44d17' 3.85"N)
    Upper Right (  347106.727,  261032.080) (116d15'59.44"W, 49d 0'18.33"N)
    Lower Right (  347106.727, -261092.845) (116d38'32.78"W, 44d19' 0.07"N)
    Center      (  -31134.301,     -30.382) (121d24'26.75"W, 46d44'56.53"N)

我在下面编写了以下探索性代码 (C++)。当我对纬度和经度执行转换时，生成的像素位置不是我所期望的。例如，当尝试将中心经度和纬度位置转换为 x、y 像素位置时。

   long = 121.407431 and  latitude = 46.749036  This is the center point from the GDALInfo output converted to decimal form.

我得到:

    x = -6234845.2270864788 and y = 5392747.8923152313  

如我所料:

    (Approx.)  x = 8935 and y = 6167

很明显，我的转换设置不正确，或者需要额外的转换和/或缩放。这是我在 C++ 中的探索性代码。这段代码来源于stack overflow的类似问题，但是结果不是我所期望的。感谢您提供的任何帮助。

    GDALAllRegister();

    GDALDataset* poDataset = (GDALDataset *)GDALOpen("D:\\programming\\gps\\Seattle SEC 94.tif", GA_ReadOnly);
    const char * szProjection = poDataset->GetProjectionRef();
    OGRSpatialReference dst(szProjection);
    OGRSpatialReference src;
    src.SetWellKnownGeogCS("WGS84");

    OGRErr error = dst.Validate();
    error = src.Validate();

    OGRCoordinateTransformation* pTransform = OGRCreateCoordinateTransformation(&src, &dst);

    // Center point of map in lattitude and longitude
    double lng = 121.407431;
    double lat = 46.749036;

    double x = lng;
    double y = lat;

    pTransform->Transform(1, &x, &y);  

x 和 x 的值不是我预期的。
实际:x = -6234845.2270864788 和 y = 5392747.8923152313
预期:(大约)x = 8935 和 y = 6167

Answer 1

需要最后一步，您必须将图像地理变换转换为像素:

double georef[6];
poDataset->GetGeoTransform(georef);

int pixelX = int((x - georef[0]) / georef[1]);
int pixelY = int((y - georef[3]) / georef[5]);