[gdal-dev] Reprojecting geometries

[Pat Blair]

For anyone else on the list who may have interest:  this seems to be an
excellent example.  I had trouble on the first try because I didn't have
PROJ.4 installed.  However, I made a quick visit to
http://www.remotesensing.org/proj/ and downloaded the libraries.  After
that, I started getting the results I was seeking.

Thank you, Mr. Warmerdam.

-----Original Message-----
From: Frank Warmerdam [mailto:warmerdam at pobox.com] 
Sent: Thursday, May 15, 2008 2:12 PM
To: Pat Blair
Cc: gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Reprojecting geometries

Pat Blair wrote:
> Hello.  I am using GDAL/OGR with Python (2.5) bindings.  I am 
> wondering if there are any ready examples that demonstrate the correct

> process for reprojecting a geometry.  I am iterating through features 
> coming from a shapefile, and the relevant code I have so far looks
similar to this...
>  
> input_ds = self.get_driver().Open(r'C:/path/to/my/shapefile.shp')
> input_layer = input_ds.GetLayer(0)
> input_feature = input_layer.GetNextFeature() while input_feature is 
> not None:
>     geometry = input_feature.GetGeometryRef()       
>     print geometry.ExportToWkt() # Looks like good WKT!  Hooray!
>     #
>     # At this point, I seem to have the geometry... but I need to
>     # reproject it to WGS-84.
>     #
>     # ...then I do some other stuff (destroying the feature, moving to

> the next one, etc.)
>  
> I have found discussion of creating spatial references in the 
> documentation and I have the idea that I can create a spatial 
> reference simply enough using the osr.SpatialReference class and 
> SetWellKnownGeogCS( "WGS84" ).  However, I'm unsure how to use it once

> I have it.  Is it possible for me to do this?

Pat,

Basically you need to create an OGRCoordinateTransformation object for
the transform from the native coordinate system to WGS84 and then call
the Transform() method on the geometry with that object.

Without testing it should be something like this:

   src_srs = input_layer.GetSpatialRef()

   dst_srs = osr.SpatialReference()
   dst_srs.SetFromUserInput( 'WGS84' )

   ct = osr.CoordinateTransformation( src_srs, dst_srs )

   ....
      geometry = input_feature.GetGeometryRef()
      print geometry.ExportToWkt()

      geometry_wgs84 = geometry.Clone()
      geometry_wgs84.Transform( ct )
      print geometry.ExportToWkt()

Keep in mind the above example creates a copy of the geometry to
transform rather than transform the one currently owned by the feature,
though I suppose you could do that too.  Just keep in mind the
Transform() method operates on the geometry it is invoked on.

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    | President OSGeo,
http://osgeo.org