[gdal-dev] Help Understanding Building Geo Transformation from scratch.

Cassanova, Bill BCassanova at weather.com
Thu Jul 30 10:21:59 EDT 2009


Hi all,

I am looking for some advice.

I have a netcdf file that contains meteorology information.  I know the
number of rows and columns of the data as well as the corner points.
The data is cylindrical-equidistant and un-projected.

I need to be able to extract data from the grid based on lat-lon
information.  Since the data is un-projected it would be easy to built a
simple ratio based on lat and lon differences between corner points and
then using simple ratio math navigate into the grid...However, I was
hoping there was a way for GDAL to do this for me.

Since the netcdf cannot be natively understood by GDAL my assumptions
are that I will need to build my own projection and create a projection
transformer.  Since I am not using GDAL to read the data source it would
seem that I can set the hSrcDS dataset to NULL.

GDALDatasetH hSrcDS, hDstDS;
char *pszDstWKT = NULL;
const char *pszSrcWKT = NULL;

hSrcDS = hDstDS = NULL;

OGRSpatialReferenceH hDEST;
hDEST = OSRNewSpatialReference(NULL);

OSRSetWellKnownGeogCS(hDEST, "WGS84");

OSRExportToWkt(hDEST, &pszDstWKT);
OSRExportToWkt(hDEST, &pszSrcWKT);

OSRDestroySpatialReference(hDEST);


void *my_transformer = GDALCreateGenImgProjTransformer(NULL,       //
Since GDAL is not reading the data-set
                                                       pszSrcWKT,  //
WGS84 is unprojected
                                                       NULL,       //
There is no output
                                                       pszDstWKT,  //
This is the big question?  What should be here?
                                                       FALSE, 0.0, 0);


// Test the corner point.  The output should be 0, 0 since this is the
origin of the grid
double x_coord = -127.0;
double y_coord = 51.0;

// Transform.

int my_result = GDALGenImgProjTransform(my_transformer, TRUE, 1,
&x_coord, &y_coord,
                              NULL, &projSuccess);



Not unexpectedly the results of x_coord, and y_coord are the original
values.  What I don't understand is what parameters should be passed
into GDALCreateGenImgProjTransformer when what I want as output is the
Cartesian coordinates in I,J format that indicate how I am supposed to
index into my 2 dimensional array of data.

Thanks in advance for the help,
Bill


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


More information about the gdal-dev mailing list