[Gdal-dev] Re: Clip shapefile into a polygon

["SAEZ Laurent - CETE Méditerr./DI/ETER"]

Frank,

Thanks for your quick answer.
I'm not a programmer...
I've wrote anyway a script in Python using Geometry.Intersection method. That was my first script in Python.
I've taken your script vec_tr.py as sample.
Miracle, my first script works good !
However, I've noted some strange things :
- Feature.SetFrom doesn't copy field values,
- Geometry.Intersection can return empty geometries : In my case, it creates plines with 0 nodes. The generated shapefile can't be open with MapInfo® 7.8.
I try to convert it with ogr2ogr in MIF/MID format. MapInfo® can't open it as well.

Thanks a lot for all. I hope that in a few months I will be able to help somebody.

Best regards.

Here is my script : test_clip.py
"
import osr
import ogr
import string
import sys
import os

#############################################################################
def Usage():
    print 'Usage: test_clip.py infile outfile clipfile'
    print
    sys.exit(1)

#############################################################################
# Argument processing.

infile = None
outfile = None
clipfile = None

i = 1
while i < len(sys.argv):
    arg = sys.argv[i]

    if infile is None:
        infile = arg

    elif outfile is None:
        outfile = arg

    elif clipfile is None:
        clipfile = arg

    else:
        Usage()

    i = i + 1

if outfile is None:
    Usage()
if clipfile is None:
    Usage()

#############################################################################
# Open the datasource to operate on.

in_ds = ogr.Open( infile, update = 0 )

in_layer = in_ds.GetLayer( 0 )

in_defn = in_layer.GetLayerDefn()

clip_ds = ogr.Open( clipfile, update = 0 )

clip_layer = clip_ds.GetLayer( 0 )
#############################################################################
#    Create output file with similar information.

shp_driver = ogr.GetDriverByName( 'ESRI Shapefile' )
if os.path.isfile(outfile) is True:
    shp_driver.DeleteDataSource( outfile )

shp_ds = shp_driver.CreateDataSource( outfile )

shp_layer = shp_ds.CreateLayer( in_defn.GetName(),
                                geom_type = in_defn.GetGeomType(),
                                srs = in_layer.GetSpatialRef() )

in_field_count = in_defn.GetFieldCount()

for fld_index in range(in_field_count):
    src_fd = in_defn.GetFieldDefn( fld_index )
    
    fd = ogr.FieldDefn( src_fd.GetName(), src_fd.GetType() )
    fd.SetWidth( src_fd.GetWidth() )
    fd.SetPrecision( src_fd.GetPrecision() )
    shp_layer.CreateField( fd )

#############################################################################
# Process all features in input layer.
# I suppose my clipfile contains only one feature !! (departemental borders)
clip_feat = clip_layer.GetNextFeature()
clip_geom = clip_feat.GetGeometryRef().Clone()
    
in_feat = in_layer.GetNextFeature()
while in_feat is not None:

    geom = in_feat.GetGeometryRef().Clone()

    geom = in_feat.GetGeometryRef().Clone().Intersection(clip_geom)
# I want to know if the geometry is empty. So, I count points;
# If the geometry is a multiline, GetPointCount return 0.
# So, I count objects.
# I know this is not rigorous
    if (geom.GetGeometryCount() + geom.GetPointCount())!=0:
        out_feat = ogr.Feature( feature_def = shp_layer.GetLayerDefn() )
        out_feat.SetFrom( in_feat )
        out_feat.SetGeometryDirectly( geom )
#=====================================================
# I copy attribute values...
        for fld_index2 in range(out_feat.GetFieldCount()):
            src_field = in_feat.GetField(fld_index2)
            out_feat.SetField(fld_index2, src_field)
#=====================================================    
        shp_layer.CreateFeature( out_feat )
        out_feat.Destroy()

    in_feat.Destroy()
    in_feat = in_layer.GetNextFeature()

#############################################################################
# Cleanup

clip_feat.Destroy()
clip_ds.Destroy()
shp_ds.Destroy()
in_ds.Destroy()

"