XML/KML - How to test if a polygon 'touches' another polygon in a bash script

Question

I am using Google My Maps to map a suburb (as a polygon) on one layer and then on a second layer, all the (rural) properties (also as polygons) that exist in and around that suburb. The polygons are using geo coordinates to define the points. Eg (somewhat redacted):

        <Polygon>
          <outerBoundaryIs>
            <LinearRing>
              <tessellate>1</tessellate>
              <coordinates>
                150.523000,-34.720000,0
                150.524000,-34.719000,0
                150.524000,-34.719000,0
                150.524000,-34.719000,0
                150.523000,-34.719000,0
                150.523000,-34.720000,0
                150.523000,-34.720000,0
              </coordinates>
            </LinearRing>
          </outerBoundaryIs>
        </Polygon>

I am then using a bash script to modify the XML/KML to meet my requirements.

I was hoping that someone would have a magic bullet that, via the command line, allows me to specify an XML file and the XPATH to two polygons (within the XML file) and have it tell me if the two polygons intersect - something like:

check_polygons "/path/to/XML" "//_:xpath/to/polygon1" "//_:xpath/to/polygon2"

and it somehow (text or exit statuses) tells me if the two polygons intersect at any point.

Any assistance is greatly appreciated.

Please read over [this guide](https://stackoverflow.com/help/asking). http://idownvotedbecau.se/noattempt/ Aside from that, I'd say try a lot of Google Fu for examples like [this](https://gis.stackexchange.com/questions/82697/shapely-linestring-and-polygon-intersect) and use `bc` or `awk` or `perl` to handle the complex math. (Better to just write the whole thing in `perl`, I suspect.) — Paul Hodges, May 29 '21 at 16:09

score 1 · Accepted Answer · answered May 30 '21 at 13:37

Thanks to Paul Hodges for giving me a kick and making me put in the effort to come up with a solution.

This is very specific to my situation which is using Google My Maps to create a map with two layers. The first layer is named 'Locality' and is a polygon that maps a suburb. The second layer is named 'Property' and contains multiple polygons that map property boundaries.

The script below takes 2 or 3 arguments. The first (optional) argument is the XML namespace (defaults to http://www.opengis.net/kml/2.2), the second argument is the name of the Placemark on the 'Property' layer and the last argument is the KML filename. The script checks if the two polygons intersect. If they do, it prints True and has a return code of 0 (zero). If they don't intersect, it prints False and has a return code of 1. A return code of 2 means the XPath couldn't be found in the XML file or there was a usage error.

#!/usr/bin/python

from lxml import etree
from shapely.geometry import Polygon
import sys


def getPoly(XmlFile, XmlNS, XPath):
    tree = etree.parse(XmlFile)
    root = tree.getroot()
    XmlNamespace = {"NS": XmlNS}
    xmlList = root.xpath(XPath, namespaces=XmlNamespace)
    if not xmlList:
        return(False)
    polyCoordList = xmlList[0].text.replace(' ', '')
    coordinates = []
    for point_text in polyCoordList.split():
      floats = point_text.split(",")
      coordinates.append((float(floats[0]), float(floats[1])))
    return(Polygon(coordinates))



if __name__ == "__main__":
    XmlNS = "http://www.opengis.net/kml/2.2"
    if len(sys.argv) == 4:
        XmlNS = sys.argv[1]
        XmlPlacemarkName = sys.argv[2]
        XmlFile = sys.argv[3]
    elif len(sys.argv) == 3:
        XmlPlacemarkName = sys.argv[1]
        XmlFile = sys.argv[2]
    else:
        print("Usage: %s [NameSpace] PlacemarkName KmlFile" % (sys.argv[0]))
        print("The default namespace is %s" % XmlNS)
        exit(2)


    XPath = "./NS:Document/NS:Folder[contains(., 'Locality')]/NS:Placemark/NS:Polygon/NS:outerBoundaryIs/NS:LinearRing/NS:coordinates"
    LocalityPoly = getPoly(XmlFile, XmlNS, XPath)
    if LocalityPoly == False:
        print("Invalid XPath - %s" % XPath)
        exit(2)
    XPath = "./NS:Document/NS:Folder[contains(., 'Property')]/NS:Placemark[contains(., '%s')]/NS:Polygon/NS:outerBoundaryIs/NS:LinearRing/NS:coordinates" % XmlPlacemarkName
    PlacemarkPoly = getPoly(XmlFile, XmlNS, XPath)
    if PlacemarkPoly == False:
        XPath = "./NS:Document/NS:Folder[contains(., 'Property')]/NS:Placemark[contains(., '%s')]/NS:MultiGeometry/NS:Polygon/NS:outerBoundaryIs/NS:LinearRing/NS:coordinates" % XmlPlacemarkName
        PlacemarkPoly = getPoly(XmlFile, XmlNS, XPath)
        if PlacemarkPoly == False:
            print("Invalid XPath - %s" % XPath)
            exit(2)
    result = LocalityPoly.intersects(PlacemarkPoly)
    print(result)
    exit(not result)

Hope that may help someone else.

Bravo. Now, if you like, we can help you tune it - but you have the satisfaction of knowing you did it on your own. :) — Paul Hodges, Jun 02 '21 at 13:10

XML/KML - How to test if a polygon 'touches' another polygon in a bash script

1 Answers1