[Gdal-dev] Problems with gtiff projections

[Roger James]

Frank

I am using a build from CVS from two days ago. I will do some further
checking. These are the results I have at the moment.


The wkt string is generated by OSR.importFromEPSG on coord code 27700

PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB
1936",DATUM["OSGB_1936",SPHEROID["Airy
1830",6377563.396,299.3249646,AUTHORITY["EPSG","7001"]],AUTHORITY["EPSG"
,"6277"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0
.0174532925199433],AUTHORITY["EPSG","4277"]],PROJECTION["Transverse_Merc
ator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-
2],PARAMETER["scale_factor",0.999601272],PARAMETER["false_easting",40000
0],PARAMETER["false_northing",-100000],UNIT["metre",1,AUTHORITY["EPSG","
9001"]],AUTHORITY["EPSG","27700"]]

The code that saves it does this

if (NULL == (pNewDataset = pDriver->Create(lpszPathName,
			m_iImagePixelWidth,
			m_iImagePixelHeight,
			iRasterCount,
			DataType,
			aOptions)))
{
	ErrorMessage.Format("Unable to create new GDAL dataset - %s",
CPLGetLastErrorMsg());
	throw (const char*)ErrorMessage;
}
#ifdef _DEBUG
DebugString.Format("SaveFile 1 - Projection %s\n", m_Wkt);
OutputDebugString(DebugString);
#endif
if (CE_None != pNewDataset->SetProjection(m_Wkt))
{
	ErrorMessage.Format("Unable to set projection - %s",
CPLGetLastErrorMsg());
	throw (const char*)ErrorMessage;
}

This is what comes back in when I re read the file.

PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB
1936",DATUM["OSGB_1936",SPHEROID["Airy
1830",6377563.396,299.3249646000044,AUTHORITY["EPSG","7001"]],AUTHORITY[
"EPSG","6277"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4277"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHO
RITY["EPSG","27700"]]

This is what listgeo says

C:\mapdata\IMAGES>listgeo test.tif
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                0                0
         0                0                0
      ModelPixelScaleTag (1,3):
         0                0                0
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,34): "OSGB 1936 / British National Grid"
      ProjectedCSTypeGeoKey (Short,1): PCS_British_National_Grid
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 27700 (name unknown)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    (      0.000,      0.000)
Lower Left    (      0.000,      0.000)
Upper Right   (      0.000,      0.000)
Lower Right   (      0.000,      0.000)
Center        (      0.000,      0.000)

This is what gdalinfo says:-

C:\mapdata\IMAGES>/vterrain/apis/gdal/bin/gdalinfo test.tif
Driver: GTiff/GeoTIFF
Size is 822, 822
Coordinate System is:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646000044,
                AUTHORITY["EPSG","7001"]],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4277"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27700"]]
Origin = (0.000000,0.000000)
Pixel Size = (0.000000,0.000000)
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Lower Left  (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Upper Right (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Lower Right (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Center      (   0.0000000,   0.0000000) (  0d 0'0.01"E,  0d 0'0.01"N)
Band 1 Block=822x3 Type=Byte, ColorInterp=Red
Band 2 Block=822x3 Type=Byte, ColorInterp=Green
Band 3 Block=822x3 Type=Byte, ColorInterp=Blue

Any ideas?

Roger

-----Original Message-----
From: gdal-dev-admin at remotesensing.org
[mailto:gdal-dev-admin at remotesensing.org] On Behalf Of Frank Warmerdam
Sent: 22 February 2003 02:02
To: gdal-dev at remotesensing.org
Subject: Re: [Gdal-dev] Problems with gtiff projections

Roger James wrote:
> Frank,
> 
>  
> 
> I create a new geotiff dataset, set the projection from wkt that has 
> projection parameters in it, and then save it. When I read it back in 
> the projection parameters are missing. For a transverse mercator these

> are things like the lat/long origin and the false easting/northing
etc.
> 
>  
> 
>  From a quick look at the code the routine GTIFGetProjTRFInfo might be
a 
> culprit as it always fails because the file "projop_wparm.csv" cannot
be 
> found.

Roger,

What version of GDAL are you using?  How are you creating and writing
the projection to the file.  I ask because in GDAL 1.1.8 there was a bug
in writing the projection to files created with the Create() method,
while
those created by the CreateCopy() were OK.

If that isn't the issue (the problem is fixed in CVS), then I will need
further details.  Generally speaking if you provide the full set of
projection parameters when writing a GeoTIFF file, then GDAL should be
able
to read and interprete the file properly when reading even if
projop_wparm.csv
isn't found.  If there is a glitch in this I would like to know, but I
will
need details of how to reproduce.  For instance, what was the WKT you
were
writing, and what is the listgeo *and* gdalinfo reports on the resulting
file?

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