Index: ogr2ogr.cpp =================================================================== RCS file: /cvs/maptools/cvsroot/gdal/ogr/ogr2ogr.cpp,v retrieving revision 1.28 diff -u -r1.28 ogr2ogr.cpp --- ogr2ogr.cpp 14 Apr 2005 14:20:24 -0000 1.28 +++ ogr2ogr.cpp 6 Jul 2005 04:13:58 -0000 @@ -117,6 +117,10 @@ #include "ogrsf_frmts.h" #include "cpl_conv.h" #include "cpl_string.h" +#include +#include +#include + CPL_CVSID("$Id: ogr2ogr.cpp,v 1.28 2005/04/14 14:20:24 fwarmerdam Exp $"); @@ -133,10 +137,16 @@ char **papszSelFields, int bAppend, int eGType ); +static void *CreateGCPTransform (const char *file); + +static OGRGeometry *warpGeometry ( + OGRGeometry *poGeometry, void *poTransform ); + static int bSkipFailures = FALSE; static int nGroupTransactions = 200; static int bPreserveFID = FALSE; static int nFIDToFetch = OGRNullFID; +static void *pTPSxform = NULL; /************************************************************************/ /* main() */ @@ -303,6 +313,10 @@ papszSelFields = CSLTokenizeStringComplex(pszSelect, " ,", FALSE, FALSE ); } + else if( EQUAL(papszArgv[iArg],"-gcp") && papszArgv[iArg+1] != NULL ) + { + pTPSxform = CreateGCPTransform(papszArgv[++iArg]); + } else if( papszArgv[iArg][0] == '-' ) { Usage(); @@ -806,6 +820,15 @@ } } + if (pTPSxform && poDstFeature->GetGeometryRef() != NULL ) { + OGRGeometry *poDstGeom = poDstFeature->StealGeometry(); + OGRGeometry *poDstGeom2; + + poDstGeom2 = warpGeometry( poDstGeom, pTPSxform ); + poDstFeature->SetGeometryDirectly(poDstGeom2); + } + + if( poDstFeature->GetGeometryRef() != NULL && bForceToPolygon ) { poDstFeature->SetGeometryDirectly( @@ -842,3 +865,113 @@ return TRUE; } +static void warpPoint (OGRPoint *poPoint, void *pXformer) +{ + double x = poPoint->getX(), + y = poPoint->getY(), + z = poPoint->getZ(); + int r = 0; + // printf("%lf %lf %lf -> ", x, y, z); + GDALTPSTransform(pXformer, FALSE, 1, &x, &y, &z, &r ); + // printf("%lf %lf %lf\n", x, y, z); + if (!r) { + printf("transform failed!\n"); + return; + } + poPoint->setX(x); + poPoint->setY(y); + poPoint->setZ(z); +} + +static void warpCurve (OGRCurve *poOut, OGRCurve *poIn, void *pXformer) +{ + OGRPoint *poPoint = new OGRPoint(); + OGRLineString *poIn2 = (OGRLineString *) poIn; + OGRLineString *poOut2 = (OGRLineString *) poOut; + for (int n = 0; n < poIn2->getNumPoints(); n++) { + poIn2->getPoint(n, poPoint); + warpPoint(poPoint, pXformer); + poOut2->addPoint(poPoint); + } + delete poPoint; +} + +static OGRGeometry *warpGeometry (OGRGeometry *poGeometry, void *poTransform) +{ + OGRGeometry *poGeomOut = NULL; + OGRwkbGeometryType eGeomType = poGeometry->getGeometryType(); + + if (eGeomType == wkbPoint) { + OGRPoint *poPoint = (OGRPoint *) poGeometry->clone(); + warpPoint(poPoint, poTransform); + poGeomOut = (OGRGeometry *)poPoint; + } + if (eGeomType == wkbLineString) { + OGRLineString *poLineIn = (OGRLineString *) poGeometry; + OGRLineString *poLineOut = new OGRLineString(); + warpCurve(poLineOut, poLineIn, poTransform); + poGeomOut = (OGRGeometry *)poLineOut; + } + else if (eGeomType == wkbPolygon) { + OGRPolygon *poPolyIn = (OGRPolygon *) poGeometry; + OGRPolygon *poPolyOut = new OGRPolygon(); + OGRLinearRing *poRingIn = poPolyIn->getExteriorRing(); + OGRLinearRing *poRingOut = new OGRLinearRing(); + + //printf("-."); + warpCurve( poRingOut, poRingIn, poTransform ); + poPolyOut->addRingDirectly( poRingOut ); + //printf(".-"); + + for (int n = 0; n < poPolyIn->getNumInteriorRings(); n++) { + poRingIn = poPolyIn->getInteriorRing(n); + poRingOut = new OGRLinearRing(); + warpCurve( poRingOut, poRingIn, poTransform ); + poPolyOut->addRingDirectly( poRingOut ); + } + + poGeomOut = (OGRGeometry *)poPolyOut; + + } + return poGeomOut; +} + +static void *CreateGCPTransform (const char *file) { + int n = 0, r = 0; + GDAL_GCP gcps[32]; + FILE *gcpfile; + void *poTransform; + + gcpfile = fopen(file, "r"); + if (!gcpfile) { + printf("Can't open GCP file %s.\n", file); + exit(0); + } + + while (!feof(gcpfile)) { + r = fscanf(gcpfile, "%lf %lf %lf %lf", + &(gcps[n].dfGCPX), &(gcps[n].dfGCPY), + &(gcps[n].dfGCPPixel), &(gcps[n].dfGCPLine) ); + if (r < 4) + continue; +/* printf( "GCP %lf %lf %lf %lf\n", + gcps[n].dfGCPX, gcps[n].dfGCPX, + gcps[n].dfGCPPixel, gcps[n].dfGCPLine ); */ + gcps[n].dfGCPZ = 0; + gcps[n].pszId = gcps[n].pszInfo = ""; + if (n++ == 32) { + printf("only 32 GCPs supported, sorry :-/\n"); + break; + } + } + fclose(gcpfile); + + poTransform = GDALCreateTPSTransformer( n, gcps, 0 ); + if (poTransform == NULL) { + printf("Couldn't generate thin plate spline transform from GCPs.\n"); + exit(1); + } + + return poTransform; +} +