xmlpolygonkmlline-intersection

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


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.


Solution

  • 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.