import math

siteboundary=None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "Site_Boundary_Auto_Theo":
        siteboundary = lyr
        break

targetlayer=None
for lyr2 in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr2.name() == "Site_Boundary_Auto_Ben":
        targetlayer = lyr2
        break

siteboundaryshapes = siteboundary.getFeatures()
targetlayershapes = targetlayer.getFeatures()

for i in siteboundaryshapes:
    siteboundarygeom = i.geometry()
    siteboundarybuffer = siteboundarygeom.buffer(1000,25)

    for j in targetlayershapes:
        targetlayergeom = j.geometry()
        targetlayeratt = j.attributes()

    if targetlayergeom.intersects(siteboundarybuffer):

        poly1 = siteboundarygeom.exportToWkt()
        poly2 = targetlayergeom.exportToWkt()

        wkt1 = []
        wkt2 = []

        poly1 = poly1.replace("Polygon ((","")
        poly2 = poly2.replace("Polygon ((","")

        poly1 = poly1.replace("Multi","")
        poly2 = poly2.replace("Multi","")

        poly1 = poly1.replace(")","")
        poly2 = poly2.replace(")","")

        poly1 = poly1.replace("(","")
        poly2 = poly2.replace("(","")

        poly1Split = poly1.split(", ")
        poly2Split = poly2.split(", ")

        for item in poly1Split:
            xysplit = item.split(" ")
            xysplit[0] = float(xysplit[0])
            xysplit[1] = float(xysplit[1])
            wkt1.append(xysplit)

        for item in poly2Split:
            xysplit = item.split(" ")
            xysplit[0] = float(xysplit[0])
            xysplit[1] = float(xysplit[1])
            wkt2.append(xysplit)

        result = []
        bresult = []

        for item in wkt1[0::1]:
            for o, j in enumerate(wkt2[0:len(wkt2)-1:1]):
                x1 = item[0]
                xy1 = wkt2[o]
                xy2 = wkt2[o+1]
                x2 = xy1[0]
                x3 = xy2[0]

                y1 = item[1]
                y2 = xy1[1]
                y3 = xy2[1]

                deltaX = (x2-x1)
                deltaY = (y2-y1)

                deltaX2 = (x3-x1)
                deltaY2 = (y3-y1)

                dist = math.sqrt((deltaX*deltaX) + (deltaY*deltaY))
                bearing1 = math.atan2(deltaY,deltaX)

                dist2 = math.sqrt((deltaX2*deltaX2) + (deltaY2*deltaY2))
                bearing2 = math.atan2(deltaY2,deltaX2)

                xmp = (x2+x3)/2
                ymp = (y2+y3)/2

                deltaXMP = (xmp-x1)
                deltaYMP = (ymp-y1)

                distMP = math.sqrt((deltaXMP*deltaXMP) + (deltaYMP*deltaYMP))

                if distMP > dist or distMP > dist2:
                    result.append(dist)
                    bresult.append(bearing1)
                else:
                    internalangle = abs(bearing2 - bearing1)
                    oppositelength = math.sqrt((((dist*dist)+(dist2*dist2))-(2*dist*dist2*math.cos(internalangle))))
                    secondinternalangle = math.asin(dist*((math.sin(internalangle))/oppositelength))
                    if bearing2 < 0:
                        polygonbearing = ((2*math.pi)+(bearing2+((math.pi/2)-secondinternalangle)))
                        bresult.append(polygonbearing)
                        polygondist = abs(dist2*math.sin(secondinternalangle))
                        result.append(polygondist)
                    else:
                        polygonbearing = bearing2-((math.pi/2)-secondinternalangle)
                        bresult.append(polygonbearing)
                        polygondist = abs(dist2*math.sin(secondinternalangle))
                        result.append(polygondist)

        for item2 in wkt2[0::1]:
            for p,q in enumerate (wkt1[0:len(wkt1)-1:1]):
                rx1 = item2[0]
                rxy1 = wkt1[p]
                rxy2 = wkt1[p+1]
                rx2 = rxy1[0]
                rx3 = rxy2[0]

                ry1 = item2[1]
                ry2 = rxy1[1]
                ry3 = rxy2[1]

                deltarX = (rx2-rx1)
                deltarY = (ry2-ry1)

                deltarX2 = (rx3-rx1)
                deltarY2 = (ry3-ry1)

                rdist = math.sqrt((deltarX*deltarX) + (deltarY*deltarY))
                rbearing1 = math.atan2(deltarY,deltarX)

                rdist2 = math.sqrt((deltarX2*deltarX2) + (deltarY2*deltarY2))
                rbearing2 = math.atan2(deltarY2,deltarX2)

                rxmp = (rx2+rx3)/2
                rymp = (ry2+ry3)/2

                deltarXMP = (rxmp-rx1)
                deltarYMP = (rymp-ry1)

                distrMP = math.sqrt((deltarXMP*deltarXMP) + (deltarYMP*deltarYMP))

                if distrMP > rdist or distrMP > rdist2:
                    result.append(rdist)
                    bresult.append(rbearing1)
                else:
                    rinternalangle = abs(rbearing2 - rbearing1)
                    roppositelength = math.sqrt((((rdist*rdist)+(rdist2*rdist2))-(2*rdist*rdist2*math.cos(rinternalangle))))
                    rsecondinternalangle = math.asin(rdist*((math.sin(rinternalangle))/roppositelength))
                    if bearing2 < 0:
                        rpolygonbearing = ((2*math.pi)+(rbearing2+((math.pi/2)-rsecondinternalangle)))
                        bresult.append(rpolygonbearing)
                        rpolygondist = abs(rdist2*math.sin(rsecondinternalangle))
                        result.append(rpolygondist)
                    else:
                        rpolygonbearing = rbearing2-((math.pi/2)-rsecondinternalangle)
                        bresult.append(rpolygonbearing)
                        rpolygondist = abs(rdist2*math.sin(rsecondinternalangle))
                        result.append(rpolygondist)

        mindist = min(result)
        bearingindex = result.index(mindist)
        minbearing = bresult[bearingindex]
        minbearing = 90 - minbearing*(180/math.pi)

        print(mindist)
        print(minbearing)
        print("Next Planning Permission")