text: #!/usr/bin/env python
  # -*- coding: utf-8 -*-
  ###############################################################################
  #
  # Project:  GDAL/OGR Test Suite
  # Purpose:  Test topological overlay methods in Layer class.
  # Author:   Ari Jolma <ari.jolma@aalto.fi>
  #
  ###############################################################################
  # Copyright (c) 2012, Ari Jolma <ari.jolma@aalto.fi>
  # Copyright (c) 2012-2013, Even Rouault <even dot rouault at mines-paris dot org>
  #
  # Permission is hereby granted, free of charge, to any person obtaining a
  # copy of this software and associated documentation files (the "Software"),
  # to deal in the Software without restriction, including without limitation
  # the rights to use, copy, modify, merge, publish, distribute, sublicense,
  # and/or sell copies of the Software, and to permit persons to whom the
  # Software is furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included
  # in all copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  # DEALINGS IN THE SOFTWARE.
  ###############################################################################
  
  import sys
  
  sys.path.append( '../pymod' )
  
  import gdaltest
  import ogrtest
  
  from osgeo import ogr
  
  ###############################################################################
  # Common usage tests.
  
  ds = None
  A = None
  B = None
  C = None
  pointInB = None
  D1 = None
  D2 = None
  empty = None
  
  def recreate_layer_C():
      global C
  
      ds.DeleteLayer( 'C' )
      C = ds.CreateLayer('C')
  
  def is_same(A,B):
  
      A.ResetReading()
      B.ResetReading()
      while True:
          fA = A.GetNextFeature()
          fB = B.GetNextFeature()
          if fA is None:
              if fB is not None:
                  return False
              return True
          if fB is None:
              return False
          if fA.Equal(fB) != 0:
              return False
  
  def algebra_setup():
  
      global ds, A, B, C, pointInB, D1, D2, empty
  
      if not ogrtest.have_geos():
          return 'skip'
  
      # Create three memory layers for intersection.
  
      ds = ogr.GetDriverByName('Memory').CreateDataSource( 'wrk' )
  
      A = ds.CreateLayer( 'A' )
      A.CreateField( ogr.FieldDefn("A", ogr.OFTInteger) )
  
      B = ds.CreateLayer( 'B' )
      B.CreateField( ogr.FieldDefn("B", ogr.OFTString) )
  
      pointInB = ds.CreateLayer( 'pointInB' )
  
      C = ds.CreateLayer( 'C' )
  
      # Add polygons.
  
      a1 = 'POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))'
      a2 = 'POLYGON((5 2, 5 3, 7 3, 7 2, 5 2))'
      b1 = 'POLYGON((2 1, 2 4, 6 4, 6 1, 2 1))'
      pointInB1 = 'POINT(3 3)'
  
      feat = ogr.Feature( A.GetLayerDefn() )
      feat.SetField('A',1)
      feat.SetGeometryDirectly( ogr.Geometry(wkt = a1) )
      A.CreateFeature( feat )
  
      feat = ogr.Feature( A.GetLayerDefn() )
      feat.SetField('A',2)
      feat.SetGeometryDirectly( ogr.Geometry(wkt = a2) )
      A.CreateFeature( feat )
  
      feat = ogr.Feature( B.GetLayerDefn() )
      feat.SetField('B','first')
      feat.SetGeometryDirectly( ogr.Geometry(wkt = b1) )
      B.CreateFeature( feat )
  
      feat = ogr.Feature( pointInB.GetLayerDefn() )
      feat.SetGeometryDirectly( ogr.Geometry(wkt = pointInB1) )
      pointInB.CreateFeature( feat )
  
      d1 = 'POLYGON((1 2, 1 3, 3 3, 3 2, 1 2))'
      d2 = 'POLYGON((3 2, 3 3, 4 3, 4 2, 3 2))'
  
      D1 = ds.CreateLayer( 'D1' )
  
      feat = ogr.Feature( D1.GetLayerDefn() )
      feat.SetGeometryDirectly( ogr.Geometry(wkt = d1) )
      D1.CreateFeature( feat )
  
      feat = ogr.Feature( D1.GetLayerDefn() )
      feat.SetGeometryDirectly( ogr.Geometry(wkt = d2) )
      D1.CreateFeature( feat )
  
      D2 = ds.CreateLayer( 'D2' )
  
      feat = ogr.Feature( D2.GetLayerDefn() )
      feat.SetGeometryDirectly( ogr.Geometry(wkt = d1) )
      D2.CreateFeature( feat )
  
      feat = ogr.Feature( D2.GetLayerDefn() )
      feat.SetGeometryDirectly( ogr.Geometry(wkt = d2) )
      D2.CreateFeature( feat )
  
      empty = ds.CreateLayer('empty')
  
      return 'success'
  
  
  def algebra_intersection():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Intersection; this should return two rectangles
  
      err = A.Intersection( B, C )
  
      if err != 0:
          gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Intersection' )
          return 'fail'
  
      C_defn = C.GetLayerDefn()
      if C_defn.GetFieldCount() != 2 \
         or C_defn.GetFieldDefn(0).GetName() != 'A' \
         or C_defn.GetFieldDefn(0).GetType() != ogr.OFTInteger \
         or C_defn.GetFieldDefn(1).GetName() != 'B' \
         or C_defn.GetFieldDefn(1).GetType() != ogr.OFTString:
          gdaltest.post_reason( 'Did not get expected output schema.' )
          return 'fail'
  
      if C.GetFeatureCount() != 2:
          gdaltest.post_reason( 'Layer.Intersection returned '+str(C.GetFeatureCount())+' features' )
          return 'fail'
  
      f1 = (ogr.Geometry(wkt = 'POLYGON ((2 3,3 3,3 2,2 2,2 3))'),1,'first')
      f2 = (ogr.Geometry(wkt = 'POLYGON ((5 2,5 3,6 3,6 2,5 2))'),2,'first')
  
      C.ResetReading();
      while 1:
          feat = C.GetNextFeature()
          if not feat: break
  
          g = feat.GetGeometryRef()
          if g.Equals(f1[0]):
              if feat.GetField('A') != f1[1] or feat.GetField('B') != f1[2]:
                  gdaltest.post_reason('Did not get expected field values.')
                  return 'fail'
          elif g.Equals(f2[0]):
              if feat.GetField('A') != f2[1] or feat.GetField('B') != f2[2]:
                  gdaltest.post_reason('Did not get expected field values.')
                  return 'fail'
          else:
              gdaltest.post_reason( 'Layer.Intersection returned wrong geometry: '+g.ExportToWkt() )
              return 'fail'
  
      # This time we test with PROMOTE_TO_MULTI and pre-created output fields.
      recreate_layer_C()
      C.CreateField( ogr.FieldDefn("A", ogr.OFTInteger) )
      C.CreateField( ogr.FieldDefn("B", ogr.OFTString) )
  
      err = A.Intersection( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if err != 0:
          gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Intersection' )
          return 'fail'
  
      if C.GetFeatureCount() != 2:
          gdaltest.post_reason( 'Layer.Intersection returned '+str(C.GetFeatureCount())+' features' )
          return 'fail'
  
      f1 = (ogr.Geometry(wkt='MULTIPOLYGON (((2 3,3 3,3 2,2 2,2 3)))'),1,'first')
      f2 = (ogr.Geometry(wkt='MULTIPOLYGON (((5 2,5 3,6 3,6 2,5 2)))'),2,'first')
  
      C.ResetReading();
      while 1:
          feat = C.GetNextFeature()
          if not feat: break
  
          g = feat.GetGeometryRef()
          if g.Equals(f1[0]):
              if feat.GetField('A') != f1[1] or feat.GetField('B') != f1[2]:
                  gdaltest.post_reason('Did not get expected field values. (1)')
                  print(feat.GetField('A'))
                  print(feat.GetField('B'))
                  return 'fail'
          elif g.Equals(f2[0]):
              if feat.GetField('A') != f2[1] or feat.GetField('B') != f2[2]:
                  gdaltest.post_reason('Did not get expected field values. (2)')
                  print(feat.GetField('A'))
                  print(feat.GetField('B'))
                  print(feat.GetField('B'))
                  return 'fail'
          else:
              gdaltest.post_reason( 'Layer.Intersection returned wrong geometry: '+g.ExportToWkt() )
              return 'fail'
  
      recreate_layer_C()
  
      # Intersection with self ; this should return 2 polygons
  
      err = D1.Intersection( D2, C )
  
      if err != 0:
          gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Intersection' )
          return 'fail'
  
      if not is_same(D1, C):
          gdaltest.post_reason( 'D1 != C' )
          return 'fail'
  
      return 'success'
  
  
  def algebra_union():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Union; this should return 5 polygons
  
      err = A.Union( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Union' )
              return 'fail'
  
          if C.GetFeatureCount() != 5:
              gdaltest.post_reason( 'Layer.Union returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.Union( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Union' )
              return 'fail'
  
          if C.GetFeatureCount() != 5:
              gdaltest.post_reason( 'Layer.Union returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # Union with self ; this should return 2 polygons
  
      err = D1.Union( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Union' )
              return 'fail'
  
          if not is_same(D1, C):
              gdaltest.post_reason( 'D1 != C' )
              return 'fail'
  
      recreate_layer_C()
  
      # Union of a polygon and a point within : should return the point and the polygon (#4772)
  
      err = B.Union( pointInB, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Union' )
              return 'fail'
  
          if C.GetFeatureCount() != 2:
              gdaltest.post_reason( 'Layer.Union returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      return 'success'
  
  
  def algebra_symdifference():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # SymDifference; this should return 3 polygons
  
      err = A.SymDifference( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.SymDifference' )
              return 'fail'
  
          if C.GetFeatureCount() != 3:
              gdaltest.post_reason( 'Layer.SymDifference returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.SymDifference( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.SymDifference' )
              return 'fail'
  
          if C.GetFeatureCount() != 3:
              gdaltest.post_reason( 'Layer.SymDifference returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # SymDifference with self ; this should return 0 features
  
      err = D1.SymDifference( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.SymDifference' )
              return 'fail'
  
      if C.GetFeatureCount() != 0:
              gdaltest.post_reason( 'Layer.SymDifference returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      return 'success'
  
  
  def algebra_identify():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Identity; this should return 4 polygons
  
      err = A.Identity( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Identity' )
              return 'fail'
  
          if C.GetFeatureCount() != 4:
              gdaltest.post_reason( 'Layer.Identity returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.Identity( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Identity' )
              return 'fail'
  
          if C.GetFeatureCount() != 4:
              gdaltest.post_reason( 'Layer.Identity returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # Identity with self ; this should return 2 polygons
  
      err = D1.Identity( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Identity' )
              return 'fail'
  
          if not is_same(D1, C):
              gdaltest.post_reason( 'D1 != C' )
              return 'fail'
  
      return 'success'
  
  def algebra_update():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Update; this should return 3 polygons
  
      err = A.Update( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Update' )
              return 'fail'
  
          if C.GetFeatureCount() != 3:
              gdaltest.post_reason( 'Layer.Update returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.Update( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Update' )
              return 'fail'
  
          if C.GetFeatureCount() != 3:
              gdaltest.post_reason( 'Layer.Update returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # Update with self ; this should return 2 polygons
  
      err = D1.Update( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Update' )
              return 'fail'
  
          if not is_same(D1, C):
              gdaltest.post_reason( 'D1 != C' )
              return 'fail'
  
      return 'success'
  
  
  def algebra_clip():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Clip; this should return 2 polygons
  
      err = A.Clip( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Clip' )
              return 'fail'
  
          if C.GetFeatureCount() != 2:
              gdaltest.post_reason( 'Layer.Clip returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.Clip( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Clip' )
              return 'fail'
  
          if C.GetFeatureCount() != 2:
              gdaltest.post_reason( 'Layer.Clip returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # Clip with self ; this should return 2 polygons
  
      err = D1.Update( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Clip' )
              return 'fail'
  
          if not is_same(D1, C):
              gdaltest.post_reason( 'D1 != C' )
              return 'fail'
  
      return 'success'
  
  
  def algebra_erase():
      if not ogrtest.have_geos():
          return 'skip'
  
      recreate_layer_C()
  
      # Erase; this should return 2 polygons
  
      err = A.Erase( B, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Erase' )
              return 'fail'
  
          if C.GetFeatureCount() != 2:
              gdaltest.post_reason( 'Layer.Erase returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      err = A.Erase( B, C, options = ['PROMOTE_TO_MULTI=YES'] )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Erase' )
              return 'fail'
  
          if C.GetFeatureCount() != 2:
              gdaltest.post_reason( 'Layer.Erase returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
      # Erase with self ; this should return 0 features
  
      err = D1.Erase( D2, C )
  
      if ogrtest.have_geos():
          if err != 0:
              gdaltest.post_reason( 'got non-zero result code '+str(err)+' from Layer.Erase' )
              return 'fail'
  
      if C.GetFeatureCount() != 0:
              gdaltest.post_reason( 'Layer.Erase returned '+str(C.GetFeatureCount())+' features' )
              return 'fail'
  
      recreate_layer_C()
  
  
      # Erase with empty layer (or no intersection)
  
      A.Erase( empty, C )
  
      if C.GetFeatureCount() != A.GetFeatureCount():
          gdaltest.post_reason( 'Layer.Erase returned '+str(C.GetFeatureCount())+' features' )
          return 'fail'
  
      A.ResetReading()
      feat_a = A.GetNextFeature()
      feat_c = C.GetNextFeature()
      if feat_a.Equal(feat_c) != 0:
          gdaltest.post_reason('features not identical')
          feat_a.DumpReadable()
          feat_c.DumpReadable()
          return 'fail'
  
      recreate_layer_C()
  
      A.Erase( empty, C, options = ['PROMOTE_TO_MULTI=YES']  )
  
      if C.GetFeatureCount() != A.GetFeatureCount():
          gdaltest.post_reason( 'Layer.Erase returned '+str(C.GetFeatureCount())+' features' )
          return 'fail'
  
      recreate_layer_C()
  
      return 'success'
  
  def algebra_cleanup():
      if not ogrtest.have_geos():
          return 'skip'
  
      global ds, A, B, C, pointInB, D1, D2, empty
  
      D2 = None
      D1 = None
      pointInB = None
      C = None
      B = None
      A = None
      empty = None
      ds = None
  
      return 'success'
  
  
  gdaltest_list = [
      algebra_setup,
      algebra_intersection,
      algebra_union,
      algebra_symdifference,
      algebra_identify,
      algebra_update,
      algebra_clip,
      algebra_erase,
      algebra_cleanup,
      ]
  
  if __name__ == '__main__':
  
      gdaltest.setup_run( 'algebra' )
  
      gdaltest.run_tests( gdaltest_list )
  
      gdaltest.summarize()
  