[gdal-dev] Question about GDALCreateGenImgProjTransformer

[Cassanova, Bill]

Hi All,

 

I am experiencing a bit of confusion regarding the use of
GDALCreateGenImgProjTransformer.

 

I am attempting to create a ProjTransformer in 2 different ways to prove
to myself that I would get the same result...The problem is I don't but
don't understand exactly why.

 

The first method I am trying is to open the GRIB2 dataset using GDALOpen
to get a GDAL dataset handle.  This handle is then passed to
GDALCreateGenImgProjTransformer as the first argument with the second
argument being passed as NULL.  The destination dataset is one created
using the well-known "WGS84" Tag...For example

 

OGRSpatialReferenceH hDEST;

hDEST = OSRNewSpatialReference(NULL);

OSRSetWellKnownGeogCS(hDEST, "WGS84");

OSRExportToWkt(hDEST, &pszDstWKT);

OSRDestroySpatialReference(hDEST);

 

void *my_transformer = GDALCreateGenImgProjTransformer(hSrcDS,

                                                          NULL,

                                                          NULL,

                                                          pszDstWKT,

                                                          FALSE, 0.0,
0);

 

The second method passes NULL for the first argument and instead passed
the WKT instead.  The destination data set is the same as created above.

 

void *my_second_transformer = GDALCreateGenImgProjTransformer( NULL,

 
GDALGetProjectionRef( hSrcDS ),

                                                                  NULL,

 
pszDstWKT,

                                                                  FALSE,

                                                                  0.0, 0
);

 

The result of GDALGetProjectionRef( hSrcDS ) for reference :

    

PROJCS["unnamed",

    GEOGCS["Coordinate System imported from GRIB file",

        DATUM["unknown",

            SPHEROID["Sphere",6371229,0]],

        PRIMEM["Greenwich",0],

        UNIT["degree",0.0174532925199433]],

    PROJECTION["Lambert_Conformal_Conic_2SP"],

    PARAMETER["standard_parallel_1",25],

    PARAMETER["standard_parallel_2",25],

    PARAMETER["latitude_of_origin",0],

    PARAMETER["central_meridian",265],

    PARAMETER["false_easting",0],

    PARAMETER["false_northing",0]]

 

Lambert Conformal Grid Information as extracted from gdalinfo:

Origin = (-4226106.996915472671390,7240925.337282469496131)

Pixel Size = (12191.000000000000000,-12191.000000000000000)

Corner Coordinates:

Upper Left  (-4226106.997, 7240925.337) (152d52'43.04"W, 54d33'55.24"N)

Lower Left  (-4226106.997, 2023177.337) (133d25'42.19"W, 12d 5'14.28"N)

Upper Right ( 3259167.003, 7240925.337) ( 49d15'27.05"W, 57d17'54.49"N)

Lower Right ( 3259167.003, 2023177.337) ( 65d 2'28.98"W, 14d12'48.02"N)

Center      ( -483469.997, 4632051.337) (100d30'19.44"W, 40d34'21.84"N)

 

 

To test both of these I pass in the known coordinates up the upper left
hand corner of the grid ( -152.877, 54.567 ) using code:

 

double x_coord = atof( argv[2] );

double y_coord = atof( argv[3] );

 

std::cout << std::fixed << "original x_coord: " << x_coord << " y_coord:
" << y_coord << std::endl;

 

result =      GDALGenImgProjTransform( my_transformer, TRUE, 1,
&x_coord, &y_coord,

                               NULL, &panSuccess);

 

  std::cout << std::fixed << "calculated x_coord: " << x_coord << "
y_coord: " << y_coord << std::endl;

 

// Reset x-coord and y-coord from command line.

  x_coord = atof( argv[2] );

  y_coord = atof( argv[3] );

 

  std::cout << std::fixed << "second original x_coord: " << x_coord << "
y_coord: " << y_coord << std::endl;

 

   esult =

      GDALGenImgProjTransform( my_second_transformer, TRUE, 1, &x_coord,
&y_coord,

                               NULL, &panSuccess);

 

   std::cout << std::fixed << "second calculated x_coord: " << x_coord
<< " y_coord: " << y_coord << std::endl;

 

The results:

 

// This result is right on the money and exactly what I expect

original x_coord: -152.877000 y_coord: 54.567000

calculated x_coord: 0.016427 y_coord: -0.011915

 

// This result is not what I expect.  I would have assumed since I
created the transformer in a different way but using the WKT information
it should have been the same.

// instead I am getting the Lambert Easting and Northing in Meters.

second original x_coord: -152.877000 y_coord: 54.567000

second calculated x_coord: -4225906.734228 y_coord: 7241070.590061

 

Am I missing a step here?

 

Thanks,

Bill

 

++++

 

William Cassanova | Senior GFS Developer | The Weather Channel |
770.226.2368 | bcassanova at weather.com

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20101213/b9852c6d/attachment.html