Some javascript/style in this post has been disabled (
why?)
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
_______________________________________________
gdal-dev mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/gdal-dev
