[Gdal-dev] Vector warping patch to ogr2ogr

Schuyler Erle schuyler at nocat.net
Wed Jul 6 01:30:16 EDT 2005


* On  5-Jul-2005 at 10:53AM PDT, Frank Warmerdam said:
> >  
> > can u help me in transformation some vector data using the polynomial Third
> > Order method. 
>
> I don't have any programs conveniently setup to apply polynomial
> transformations to vector data; however, I do know of at least one
> person who has successfully modified the ogr2ogr code to accomplish
> this.  It is Schuyler Erle (cced on this message).  You might want to 
> ask him for a drop of his modified code. 

Hi guys, attached is a patch to ogr2ogr that allows you to specify
GCPs via standard input using the -gcp flag, which will then use the
thin plate spline warper from GDAL to warp the output points into the
matching projection. You'll probably want to use the -t_srs flag as
well. The format of the GCPs should be "x y pixel line", whitespace
separated, one GCP per line. If you need polynomial warping, you'll
want to hack the patched version. I don't think this patch is really
fit to be applied to the main distribution version of ogr2ogr, but...

Really, what I'd love to see happen is some kind of "filter"
interface, which would allow us to add an additional processing step
via command-line options, which would call out to a .so plugin, or
moral equivalent. Different vector warping modules and line
simplification are just two possible uses that come to mind. Of course
I can already imagine Frank breaking out in cold sweats reading this
email... :-) But I'd thought I'd open the can of worms and see where
they lead...

SDE


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 <gdal_alg.h>
+#include <gdal.h>
+#include <stdio.h>
+
 
 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;
+}
+



More information about the Gdal-dev mailing list