[gdal-dev] Reprojecting geometries

[Frank Warmerdam]

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