[gdal-dev] General Questions about gdal usage with Grib1 Data Set

[Cassanova, Bill]

Thank you for your help.

Couple more questions:

(1)

After the call to  GDALGenImgProjTransform( ..., ..., ...,  double *x,
double *y, ..., ... );

I am assuming that the x and y arguments contain the new lat-lon values
from the new projection and does not contain the data from the grid.  In
other words I would need to maintain the data in a matrix of my own and
then given the results of X and Y determine an appropriate offset into
the rows and columns of that matrix.  That being said, I would assume
that when data is re-projected that the raster now has a new set of
corner points but the same number of rows and columns of data as was
original.  Am I correct? Do I have anyway of determining the new corner
points?  Has the original number of rows and columns changed?

(2)  My assumption is the C API GDALCreateGenImgProjTransformer and
GDALGenImgProjTransform function calls do not have C++ equivalents
buried in one of the classes somewhere.  I looked briefly at gdalwarper
and saw a few methods that looked similar but nothing that matches the
functionality of these calls.


Thanks,
Bill

-----Original Message-----
From: Frank Warmerdam [mailto:warmerdam at pobox.com] 
Sent: Monday, July 20, 2009 3:52 PM
To: Cassanova, Bill
Cc: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] General Questions about gdal usage with Grib1
Data Set

Cassanova, Bill wrote:
> Hi,
> 
> I am new to GDal and just beginning to dig into the intricacies of the

> library and would like some guidance on the best way to proceed.
> 
> I have a data set that is in a grib1 format but with a Lambert
Conformal 
> Projection.  If need be this file can be converted to either grib2 or 
> shape file using an application called degrib.
> 
> My task is given this data-set I need to be able to extract values
from 
> the file for a given lat-lon and don't necessarily want to code the 
> intricacies of the Lambert-Conformal-SP2 projection.  So in other
words 
> given a lat-lon of 38.59, - 97.61 I want to be able to interrogate the

> raster and extract a value for that given coordinate.
> 
> There seem to be a few options and just need advice from the season
GDAL 
> guru's of the best way to proceed:
> 
> (1)     Convert the lambert-conformal grib1 to a cylindrical
equidistant 
> grib1 or tiff using gdalwarp and do the simpler ratio based math on
the 
> rectangular grid as opposed to the projected grid.
> 
> (2)     Use OGRCreateCoordinateTransformation specifying the source as

> Lambert Conformal and the destination as cylindrical equidistant.  My 
> presumption at this point is that the 
> OGRCoordinateTransformation->Transform( 1, &x, &y ) method call could 
> take a longitude as the x-argument, and latitude as the y-argument.  A

> successful transform call would replace the original lat-lon with the 
> new lat-lon and from this point the data grid ( I, J ) offsets could
be 
> calculated using simply ratios.
> 
> (3)     Are they any more direct simpler ways of accomplishing the
same 
> thing? 

Bill,

I would suggest a variation on (2).  You can call
GDALCreateGenImgProjTransformer to create a transformer from image
pixel/line
to long/lat.  The call is:

void CPL_DLL *
GDALCreateGenImgProjTransformer( GDALDatasetH hSrcDS, const char
*pszSrcWKT,
                                  GDALDatasetH hDstDS, const char
*pszDstWKT,
                                  int bGCPUseOK, double
dfGCPErrorThreshold,
                                  int nOrder );

But you can leave the hDstDS as NULL, and provide the WGS84 pszDstWKT
yourself.  If the source file coordinate system is already known by GDAL
(check gdalinfo report) then you could leave pszSrcWKT NULL and it will
be picked up from the dataset.  The bGCPUseOK can be FALSE, and the
GCPErrorThreshold 0.0, and nOrder 0.

The pointer you get back can be used with GDALGenImgProjTransform():

int CPL_DLL GDALGenImgProjTransform(
     void *pTransformArg, int bDstToSrc, int nPointCount,
     double *x, double *y, double *z, int *panSuccess );

In your case you want to pass in long/lat, and get back image
pixel/line so you would pass bDstToSrc=TRUE.

This approach is somewhat easier than OGRCreateCoordinateTransformation
because the GenImgProj stuff knows also how to apply the geotransform
step.

Best regards,
-- 
---------------------------------------+--------------------------------
------
I set the clouds in motion - turn up   | Frank Warmerdam,
warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent