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 <norm...@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.
smime.p7s
Description: S/MIME Cryptographic Signature
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
