Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm a GIS n00b. I have some LinearRings in a KML file that give the areas for a bunch of different regions. And then I (will) have some lat lons that I want to test to see if they're in any of the (about 8) different regions that I have LinearRings for. I want to use python and shapely looks like it's a good fit, but it doesn't deal with kml data at all. Is there a standard way to do things like load KML that I can then import into shapely? I'm having trouble putting all the pieces together.

share|improve this question
add comment

1 Answer

up vote 1 down vote accepted

Both Keytree and FastKML parse KML placemarks into GeoJSON-ish objects that are easy to use with Shapely. Here's an example of using Keytree:

from urllib import urlopen
from xml.etree import ElementTree

import keytree
from shapely.geometry import Point, shape

# Parse the KML doc
doc = urlopen("http://pleiades.stoa.org/places/638753/kml").read()
tree = ElementTree.fromstring(doc)
kmlns = tree.tag.split('}')[0][1:]

# Find all Polygon elements anywhere in the doc
elems = tree.findall(".//{%s}Polygon" % kmlns)

# Here's our point of interest
p = Point(28.722144580890763, 37.707799701548467)

# Filter polygon elements using this lambda (anonymous function)
# keytree.geometry() makes a GeoJSON-like geometry object from an
# element and shape() makes a Shapely object of that.
hits = filter(
    lambda e: shape(keytree.geometry(e)).contains(p),
    elems )

print hits
# Output:  [<Element {http://www.opengis.net/kml/2.2}Polygon at ...>]

Once you've found the elements that contain the point, you can use the xml module's etree API to inspect them further.

share|improve this answer
    
I was wrestling with this all afternoon today, but then I finally figured out that keytree.geometry() is pulling things in as lon,lat not lat,lon. when I switched my code so that I was encoding points like so: p = Point(lon, lat) then it finally worked for me. Thanks, so much, Sean! –  Paul Huff Sep 9 '12 at 2:31
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.