[gdal-dev] geotiff projection not showing up
Norman Goldstein
normvcr at telus.net
Mon Nov 11 11:46:47 PST 2013
Frank and Trent,
I have created a proper GeoTIFF file with the code
listed, below. There is still one oddity:
The code sets the line/sample --> x/y affine transformation
using an array of 6 doubles. The code contains a conditional
compilation that chooses between two different arrays of
6 doubles.
The #if 1 array produces a full listgeo dump.
The #else listgeo dump is missing the corner coordinates.
The difference between the two arrays is
-- The false northing for #if 1 is 0.001 .
-- The false northing for #else is 0
Is this a quirk of GDAL? of listeo? Or, should map x,y
always be set up to take on only positive values?
/////////////////////// full source code ///////////////////////////
#include <iostream>
#include "ogr_spatialref.h"
#include "gdal_priv.h"
int main( int argc, char* argv[] )
{
using namespace std;
GDALAllRegister();
GDALDriverManager* dmgr = GetGDALDriverManager();
GDALDriver* drv =
dmgr->GetDriverByName( "GTiff" );
if( nullptr == drv )
{
cerr << "Error from GetDriverByName" << endl;
return -3;
}
GDALDataset* dataset = drv->Create( "dem.tif",
300,
200,
1, // Num bands
GDT_Float32,
nullptr );
if( nullptr == dataset )
{
cerr << "Error from Create" << endl;
return -4;
}
OGRSpatialReference oSRS;
oSRS.SetProjCS( "NoWhere" );
oSRS.SetWellKnownGeogCS( "WGS84" );
oSRS.SetEquirectangular( 0.0, // Centre lat
0.0, // Centre lon
0.0, // False Easting
0.0 ); // False Northing
char* wkt = nullptr;
if( OGRERR_NONE != oSRS.exportToWkt( &wkt ) )
{
cerr << "Error from exportToWkt" << endl;
return -5;
}
cout << "WKT= " << wkt << endl;
#if 1
double coeffs[] = { 0.0, 1.0, 0.0,
0.001, 0.0, -1.0 };
#else
double coeffs[] = { 0.0, 1.0, 0.0,
0.0, 0.0, -1.0 };
#endif
if( CE_None != dataset->SetGeoTransform( coeffs ) )
{
cerr << "Error from SetGeoTransform" << endl;
return -6;
}
if( CE_None != dataset->SetProjection( wkt ) )
{
cerr << "Error from SetProjection" << endl;
return -7;
}
OGRFree( wkt );
wkt = nullptr;
GDALClose( dataset );
dataset = nullptr;
}// main
///////////////////////////////////////////////////////
############# full listgeo dump ###############
Geotiff_Information:
Version: 1
Key_Revision: 1.0
Tagged_Information:
ModelTiepointTag (2,3):
0 0 0
0 0.001 0
ModelPixelScaleTag (1,3):
1 1 0
End_Of_Tags.
Keyed_Information:
GTModelTypeGeoKey (Short,1): ModelTypeProjected
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
GTCitationGeoKey (Ascii,8): "NoWhere"
GeographicTypeGeoKey (Short,1): GCS_WGS_84
GeogCitationGeoKey (Ascii,7): "WGS 84"
GeogAngularUnitsGeoKey (Short,1): Angular_Degree
GeogSemiMajorAxisGeoKey (Double,1): 6378137
GeogInvFlatteningGeoKey (Double,1): 298.257224
ProjectedCSTypeGeoKey (Short,1): User-Defined
ProjectionGeoKey (Short,1): User-Defined
ProjCoordTransGeoKey (Short,1): CT_Equirectangular
ProjLinearUnitsGeoKey (Short,1): Linear_Meter
ProjStdParallel1GeoKey (Double,1): 0
ProjFalseEastingGeoKey (Double,1): 0
ProjFalseNorthingGeoKey (Double,1): 0
ProjCenterLongGeoKey (Double,1): 0
ProjCenterLatGeoKey (Double,1): 0
End_Of_Keys.
End_Of_Geotiff.
Projection Method: CT_Equirectangular
ProjCenterLatGeoKey: 0.000000 ( 0d 0' 0.00"N)
ProjCenterLongGeoKey: 0.000000 ( 0d 0' 0.00"E)
ProjFalseEastingGeoKey: 0.000000 m
ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/ 0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)
Corner Coordinates:
Upper Left ( 0.000, 0.001) ( 0d 0' 0.00"E, 0d 0' 0.00"N)
Lower Left ( 0.000, -199.999) ( 0d 0' 0.00"E, 0d 0' 6.47"S)
Upper Right ( 300.000, 0.001) ( 0d 0' 9.70"E, 0d 0' 0.00"N)
Lower Right ( 300.000, -199.999) ( 0d 0' 9.70"E, 0d 0' 6.47"S)
Center ( 150.000, -99.999) ( 0d 0' 4.85"E, 0d 0' 3.23"S)
########################################
On 11/09/2013 03:14 PM, Trent Piepho wrote:
> On Sat, Nov 9, 2013 at 1:07 PM, Norman Goldstein <normvcr at telus.net> wrote:
>> Things are better, now, but not quite there for me.
>> Still not able to transform pixel/line to PCS space.
>> (the listgeo dump is, below)
>>
>> I think the problem is that there is no definition of the
>>
>> ModelPixelScaleTag
>>
>> It seems that this tag, together with the ModelTiepointTag,
>> is how an affine transformation is inferred. Or by directly
>> setting
>>
>> ModelTransformationTag
>>
>> which I could do with GDALDataset's SetGeoTransform()
>> method (for defining 2D affine transformations).
> I've found that unless you call SetGeoTransform() and give an affine
> transform, most apps, including listgeo and gdalsrsinfo, aren't
> entirely happy with the georeferencing. GDAL has a function that
> computes a transform from GCPs, but it needs to be part of the GDAL
> code for the dataset driver. GDAL doesn't automatically do it when a
> user of GDAL wants an affine transform from a dataset.
>
> It seems like most code that tries to find corner coordinates and/or
> the pixel size of a raster expects an affine transform plus a
> projection. If there are a set of GCPs, then that code won't work.
> The information may well be there, like your three GCP points, but the
> GDAL user needs totally different code make use of it. It's the same
> reason most of NOAA's nautical charts don't work with apps that use
> GDAL or GDAL created GeoTIFFs. I've just started looking at the GDAL
> code, but it doesn't seem like GDAL abstracts this enough to have
> broad compatibility.
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 4243 bytes
Desc: S/MIME Cryptographic Signature
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20131111/88391feb/attachment.bin>
More information about the gdal-dev
mailing list