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