#!/usr/bin/env pytest
# -*- 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 ogrtest

from osgeo import ogr
import pytest

###############################################################################
# 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 print_layer(A):
    A.ResetReading()
    while True:
        f = A.GetNextFeature()
        if f is None:
            return
        g = f.GetGeometryRef()
        print(g.ExportToWkt())


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 test_algebra_setup():

    global ds, A, B, C, pointInB, D1, D2, empty

    if not ogrtest.have_geos():
        pytest.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')


def test_algebra_intersection():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Intersection; this should return two rectangles

    err = A.Intersection(B, C)

    assert err == 0, \
        ('got non-zero result code ' + str(err) + ' from Layer.Intersection')

    C_defn = C.GetLayerDefn()
    assert C_defn.GetFieldCount() == 2 and C_defn.GetFieldDefn(0).GetName() == 'A' and C_defn.GetFieldDefn(0).GetType() == ogr.OFTInteger and C_defn.GetFieldDefn(1).GetName() == 'B' and C_defn.GetFieldDefn(1).GetType() == ogr.OFTString, \
        'Did not get expected output schema.'

    assert C.GetFeatureCount() == 2, \
        ('Layer.Intersection returned ' + str(C.GetFeatureCount()) + ' features')

    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]):
            assert feat.GetField('A') == f1[1] and feat.GetField('B') == f1[2], \
                'Did not get expected field values.'
        elif g.Equals(f2[0]):
            assert feat.GetField('A') == f2[1] and feat.GetField('B') == f2[2], \
                'Did not get expected field values.'
        else:
            pytest.fail('Layer.Intersection returned wrong geometry: ' + g.ExportToWkt())

    # 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'])

    assert err == 0, \
        ('got non-zero result code ' + str(err) + ' from Layer.Intersection')

    assert C.GetFeatureCount() == 2, \
        ('Layer.Intersection returned ' + str(C.GetFeatureCount()) + ' features')

    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]):
            assert feat.GetField('A') == f1[1] and feat.GetField('B') == f1[2], \
                'Did not get expected field values. (1)'
        elif g.Equals(f2[0]):
            assert feat.GetField('A') == f2[1] and feat.GetField('B') == f2[2], \
                'Did not get expected field values. (2)'
        else:
            pytest.fail('Layer.Intersection returned wrong geometry: ' + g.ExportToWkt())

    recreate_layer_C()

    # Intersection with self ; this should return 2 polygons

    err = D1.Intersection(D2, C, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])

    assert err == 0, \
        ('got non-zero result code ' + str(err) + ' from Layer.Intersection')

    assert is_same(D1, C), 'D1 != C'


def test_algebra_KEEP_LOWER_DIMENSION_GEOMETRIES():
    if not ogrtest.have_geos():
        pytest.skip()

    driver = ogr.GetDriverByName('MEMORY')
    ds = driver.CreateDataSource('ds')
    layer1 = ds.CreateLayer('layer1')
    layer2 = ds.CreateLayer('layer2')

    g1 = "POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))"
    geom1 = ogr.CreateGeometryFromWkt(g1)
    feat1 = ogr.Feature(layer1.GetLayerDefn())
    feat1.SetGeometry(geom1)
    layer1.CreateFeature(feat1)

    g2 = "POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))"
    geom2 = ogr.CreateGeometryFromWkt(g2)
    feat2 = ogr.Feature(layer2.GetLayerDefn())
    feat2.SetGeometry(geom2)
    layer2.CreateFeature(feat2)

    g1 = "LINESTRING (0 0, 1 0)"
    geom1 = ogr.CreateGeometryFromWkt(g1)
    feat1 = ogr.Feature(layer1.GetLayerDefn())
    feat1.SetGeometry(geom1)
    layer1.CreateFeature(feat1)

    g2 = "LINESTRING (1 0, 2 0)"
    geom2 = ogr.CreateGeometryFromWkt(g2)
    feat2 = ogr.Feature(layer2.GetLayerDefn())
    feat2.SetGeometry(geom2)
    layer2.CreateFeature(feat2)

    layer3 = ds.CreateLayer('layer3a')
    layer1.Intersection(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])
    assert layer3.GetFeatureCount() == 0, \
        'Lower dimension geometries not removed in intersection'

    layer3 = ds.CreateLayer('layer3b')
    layer1.Intersection(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=YES'])
    assert layer3.GetFeatureCount() == 2, \
        'Lower dimension geometries not kept in intersection'

    layer3 = ds.CreateLayer('layer3c')
    layer1.Union(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])
    assert layer3.GetFeatureCount() == 4, \
        'Lower dimension geometries not removed in union'

    layer3 = ds.CreateLayer('layer3d')
    layer1.Union(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=YES'])
    assert layer3.GetFeatureCount() == 6, 'Lower dimension geometries not kept in union'

    layer3 = ds.CreateLayer('layer3e')
    layer1.Identity(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])
    assert layer3.GetFeatureCount() == 2, \
        'Lower dimension geometries not removed in identity'

    layer3 = ds.CreateLayer('layer3f')
    layer1.Identity(layer2, layer3, ['KEEP_LOWER_DIMENSION_GEOMETRIES=YES'])
    assert layer3.GetFeatureCount() == 4, \
        'Lower dimension geometries not kept in identity'


def test_algebra_union():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Union; this should return 5 polygons

    err = A.Union(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Union')

        assert C.GetFeatureCount() == 5, \
            ('Layer.Union returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.Union(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Union')

        assert C.GetFeatureCount() == 5, \
            ('Layer.Union returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Union with self ; this should return 2 polygons

    err = D1.Union(D2, C, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Union')

        assert is_same(D1, C), 'D1 != C'

    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():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Union')

        assert C.GetFeatureCount() == 2, \
            ('Layer.Union returned ' + str(C.GetFeatureCount()) + ' features')

    

def test_algebra_symdifference():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # SymDifference; this should return 3 polygons

    err = A.SymDifference(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.SymDifference')

        assert C.GetFeatureCount() == 3, \
            ('Layer.SymDifference returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.SymDifference(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.SymDifference')

        assert C.GetFeatureCount() == 3, \
            ('Layer.SymDifference returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # SymDifference with self ; this should return 0 features

    err = D1.SymDifference(D2, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.SymDifference')

    assert C.GetFeatureCount() == 0, \
        ('Layer.SymDifference returned ' + str(C.GetFeatureCount()) + ' features')


def test_algebra_identify():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Identity; this should return 4 polygons

    err = A.Identity(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Identity')

        assert C.GetFeatureCount() == 4, \
            ('Layer.Identity returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.Identity(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Identity')

        assert C.GetFeatureCount() == 4, \
            ('Layer.Identity returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Identity with self ; this should return 2 polygons

    err = D1.Identity(D2, C, ['KEEP_LOWER_DIMENSION_GEOMETRIES=NO'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Identity')

        assert is_same(D1, C), 'D1 != C'

    

def test_algebra_update():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Update; this should return 3 polygons

    err = A.Update(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Update')

        assert C.GetFeatureCount() == 3, \
            ('Layer.Update returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.Update(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Update')

        assert C.GetFeatureCount() == 3, \
            ('Layer.Update returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Update with self ; this should return 2 polygons

    err = D1.Update(D2, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Update')

        assert is_same(D1, C), 'D1 != C'

    

def test_algebra_clip():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Clip; this should return 2 polygons

    err = A.Clip(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Clip')

        assert C.GetFeatureCount() == 2, \
            ('Layer.Clip returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.Clip(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Clip')

        assert C.GetFeatureCount() == 2, \
            ('Layer.Clip returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Clip with self ; this should return 2 polygons

    err = D1.Update(D2, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Clip')

        assert is_same(D1, C), 'D1 != C'

    

def test_algebra_erase():
    if not ogrtest.have_geos():
        pytest.skip()

    recreate_layer_C()

    # Erase; this should return 2 polygons

    err = A.Erase(B, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Erase')

        assert C.GetFeatureCount() == 2, \
            ('Layer.Erase returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    err = A.Erase(B, C, options=['PROMOTE_TO_MULTI=YES'])

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Erase')

        assert C.GetFeatureCount() == 2, \
            ('Layer.Erase returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Erase with self ; this should return 0 features

    err = D1.Erase(D2, C)

    if ogrtest.have_geos():
        assert err == 0, ('got non-zero result code ' + str(err) + ' from Layer.Erase')

    assert C.GetFeatureCount() == 0, \
        ('Layer.Erase returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()

    # Erase with empty layer (or no intersection)

    A.Erase(empty, C)

    assert C.GetFeatureCount() == A.GetFeatureCount(), \
        ('Layer.Erase returned ' + str(C.GetFeatureCount()) + ' features')

    A.ResetReading()
    feat_a = A.GetNextFeature()
    feat_c = C.GetNextFeature()
    if feat_a.Equal(feat_c) != 0:
        feat_a.DumpReadable()
        feat_c.DumpReadable()
        pytest.fail('features not identical')

    recreate_layer_C()

    A.Erase(empty, C, options=['PROMOTE_TO_MULTI=YES'])

    assert C.GetFeatureCount() == A.GetFeatureCount(), \
        ('Layer.Erase returned ' + str(C.GetFeatureCount()) + ' features')

    recreate_layer_C()


def test_algebra_cleanup():
    if not ogrtest.have_geos():
        pytest.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



