[Gdal-dev] OGR Python examples

[Frank Warmerdam]

Tyler Mitchell wrote:
> 
> 
> 
> I'm starting into OGR with python and am having fun reading ogr.py and
> trying out my own little examples.
> I've got most of what I'd call the basics down pat - opening data source,
> reading a layer, feature and geometries.
> 
> Now I'm ready to write geometries and am missing something.  I've got a two
> pairs of coordinates representing a rectangle and I want to create a
> polygon out of it and save it to a file.  When I READ geometries I can
> figure out that you open a data source->layer->feature->geometry.  But I
> can't figure out what objects I need to create first to write to.  I've got
> a data source and layer created, but what next?
> 
> Can anyone give me a small example?
> 
> I'm trying to stick as close to ogr.py as possible without mixing up too
> much with other pymods, but I will consider them.
> 
> I'll post my little tutorial afterwards.

Tyler,

There are some examples in gdal/pymod/samples/get_soundg.py, and
gdal/pymod/samples/tigerpoly.py.  Also, the GDAL test suite has many
good examples:

   http://www.remotesensing.org/cgi-bin/cvsweb.cgi/osrs/gdalautotest/

In particular, the gdalautotest/ogr/ogr_oci.py is a decent example of
writing features to a layer after reading from a file, and making them up
from scratch.

This code for (from ogr_oci_4()) creates a feature in memory, and then
set's its geometry from WKT read from a file, and writes it to the existing
Oracle layer (gdaltest.oci_lyr).  The same features is written many time with
different geometries:

     dst_feat = ogr.Feature( feature_def = gdaltest.oci_lyr.GetLayerDefn() )
     wkt_list = [ '10', '2', '1', '3d_1', '4', '5', '6' ]

     for item in wkt_list:

         wkt = open( 'data/wkb_wkt/'+item+'.wkt' ).read()
         geom = ogr.CreateGeometryFromWkt( wkt )

         ######################################################################
         # Write geometry as a new Oracle feature.

         dst_feat.SetGeometryDirectly( geom )
         dst_feat.SetField( 'PRFEDEA', item )
         gdaltest.oci_lyr.CreateFeature( dst_feat )

         ######################################################################
         # Read back the feature and get the geometry.

         gdaltest.oci_lyr.SetAttributeFilter( "PRFEDEA = '%s'" % item )
         feat_read = gdaltest.oci_lyr.GetNextFeature()
         geom_read = feat_read.GetGeometryRef()

         if ogrtest.check_feature_geometry( feat_read, geom ) != 0:
             return 'fail'

         feat_read.Destroy()

     dst_feat.Destroy()

So, the key is to create a features with the ogr.Feature(<Featuredefn>) call.
Then set any attributes or geometry you might want on the feature using
SetField() and SetGeometry() (or SetGeometryDirectly() if you want the feature
to take ownership of an existing geometry).

I hope this helps.

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    | Geospatial Programmer for Rent