[gdal-dev] Building a Mercator Image from Un-projected data.. Major confusion and help requested.

[Cassanova, Bill]

Hi all,

I am looking to gain a better understanding of how I can use GDAL to
build a Mercator Image from un-projected data.

The data I have is contained in proprietary format.  For all intensive
purposes it is a grid consisting of 2464 rows by 5796 columns of data.
The origin of the data is in the upper left hand corner of the grid with
grid lat-lon dimensions of -125.99, 49.7186 for the upper left and
-65.656, 24.0938 for the lower right.

I have tried several things but not sure I get the connection well
enough to complete the project. So...Here is my sample program with
Excerpts for comments. 

std::pair< double, double > upper_left = std::make_pair( -125.99,
49.7186 );
std::pair< double, double > upper_right = std::make_pair( -65.656,
49.7186 );
std::pair< double, double > lower_left = std::make_pair( -125.99,
24.0938 );
std::pair< double, double > lower_right = std::make_pair( -65.656,
49.7186 );

int my_points = 1;
int result = 0;

double x, y;

OGRSpatialReferenceH mercator_proj = OSRNewSpatialReference( "" );
OSRSetProjCS( mercator_proj, "Mercator Projection" );
OSRImportFromProj4( mercator_proj, "+proj=merc" );

OGRSpatialReferenceH wgs84_proj = OSRNewSpatialReference( "" );
OSRSetWellKnownGeogCS( wgs84_proj,  "WGS84" );


OGRCoordinateTransformationH my_transformation =
OCTNewCoordinateTransformation( wgs84_proj, mercator_proj );

x = upper_left.first;
y = upper_left.second;

OCTTransform( my_transformation, my_points, &x, &y, NULL );

std::cout << std::fixed;

std::cout << result << " " << upper_left.first << ", " <<
upper_left.second << " = " << x << ", " << y << std::endl;

The output for my example is:
0 -125.990000, 49.718600 = -14025142.645045, 6365068.664114

The values output in meters and are exactly what I would expect and are
similar to proj4 output.  
My question is how do I convert this into actual i-j coordinates for the
image...And even more so how do I convert these to i-j coordinates
For the data since it is un-projected?

Another Method is using GDALCreateGenImgProfTransformer:

char *pszDstWKT = NULL;
char *pszSrcWKT = NULL;

OSRExportToWkt( mercator_proj, &pszDstWKT );
OSRExportToWkt( wgs84_proj, &pszSrcWKT );


void *hTransformArg;

hTransformArg = GDALCreateGenImgProjTransformer( NULL, pszSrcWKT, NULL,
pszDstWKT, FALSE, 0, 1 );


x = upper_left.first;
y = upper_left.second;

int success;

result = GDALGenImgProjTransform(hTransformArg, TRUE, 1, &x, &y,
                              NULL, &success);


std::cout << std::fixed;

td::cout << result << " " << upper_left.first << ", " <<
upper_left.second << " = " << x << ", " << y << std::endl;
// 1 -125.990000, 49.718600 = -0.001132, 0.000450
// OK, this seems sensible but is this the i-j offset into the data or
the image or both

x = lower_right.first;
y = lower_right.second;


result = GDALGenImgProjTransform(hTransformArg, TRUE, 1, &x, &y,
                              NULL, &success);


std::cout << std::fixed;
std::cout << result << " " << lower_right.first << ", " <<
lower_right.second << " = " << x << ", " << y << std::endl;
// 1 -65.656000, 49.718600 = -0.000590, 0.000450
// This is NOT Right so apparently I am doing something wrong or not
understanding the concept here

Thanks,
Bill

 

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