[gdal-dev] (no subject)

[Frank Warmerdam]

Sjur Kolberg wrote:
>  
> Hello, GDAL list;
>  
> Can anyone suggest what prevents the following GDALSuggestedWarpOutput 
> from calculating an approximate geotransform? The resulting suggested 
> GeoTransform is just a copy of the original GeoTransform.
>  
> Version is GDAL 1.6.0beta1 under Windows XP, I use Visual Studio 2005, 
> this is a MFC application. poDataset is a 9900 band GRIB 1 data set 
> from ECMWF. It appears to be read in OK, I can extract 
> plausible metadata both from the dataset and from each band. 
>  
>    SourceWKT is (CString)
> GEOGCS["Coordinate System imported from GRIB file",DATUM["unknown",
> SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
>    SourceWKT is from poDataset->GetProjectionRef()
>  
>    TargetWKT is (CString)
> PROJCS["UTM 32 North (WGS84)",GEOGCS["WGS 84",DATUM["WGS_1984",
> SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
> TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",
> 0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,
> AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],
> AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],
> PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],
> PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],
> PARAMETER["false_northing",0],UNIT["Meter",1]]
>    TargetWKT is from 
> ogrspatref.SetProjCS(msg);                        // CString msg is just a name
> ogrspatref.SetWellKnownGeogCS( "WGS84" );
> if (ogrspatref.SetUTM( zone, i ) != OGRERR_NONE)  // zone==32, i==1 for northern hem.
>  
>  
> And the code goes:
> 
> 
> void *hTransformArg = GDALCreateGenImgProjTransformer( 
> poDataset, sourceWKT, NULL, targetWKT, FALSE, 0, 1 );
> 
> CPLAssert( hTransformArg != NULL );
> 
> if (GDALSuggestedWarpOutput(pSrcDS, GDALGenImgProjTransform, 
> hTransformArg, SugGeoTransform, &nPixels, &nLines )  != CE_None )

Sjur,

I tried this (test.wkt is the UTM WKT) and it seems to work plausibly:

testepsg -t 'GEOGCS["Coordinate System imported from GRIB 
file",DATUM["unknown",SPHEROID["Sphere",6367470,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]' 
test.wkt 4 65
OGRCT: Source: +proj=longlat +a=6367470 +b=6367470 +no_defs
OGRCT: Target: +proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
(4.000000,65.000000,0.000000) -> (264409.202311,7217779.091236,0.000000)

So, the reprojection operation can apparently work.  The behavior you
are seeing appears to be congruent with reprojection doing nothing.
Is there any possibility that the PROJ library is not built into your
GDAL?

I see no other reason things should not work as you have constructed
them.  All I can see is to either debug through and see what is going
wrong, or to submit a ticket with info on how to reproduce the problem.

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