Skip to content

Instantly share code, notes, and snippets.

@misja
Created March 16, 2014 11:28
Show Gist options
  • Save misja/9581848 to your computer and use it in GitHub Desktop.
Save misja/9581848 to your computer and use it in GitHub Desktop.
Get a country for a longitude / latitude point. World Borders Dataset to be found at http://thematicmapping.org/downloads/world_borders.php
import rtree
import cPickle
import osgeo.ogr
import shapely.wkb
import shapely.speedups
shapely.speedups.enable()
class FastRtree(rtree.Rtree):
def dumps(self, obj):
return cPickle.dumps(obj, -1)
def loader():
fname = './data/TM_WORLD_BORDERS-0.3.shp'
shapefile = osgeo.ogr.Open(fname)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
geometry = feature.GetGeometryRef()
shape = shapely.wkb.loads(geometry.ExportToWkb())
obj = {
'code': feature.GetField('ISO2').lower(),
'name': feature.GetField('NAME'),
'shape': shape
}
yield (i, shape.bounds, obj)
index = FastRtree(loader())
def detect(point):
for hit in index.intersection((point.x, point.y), objects="raw"):
if hit['shape'].contains(point):
return hit
if __name__ == '__main__':
points = [
# lon / lat pairs
('nl', shapely.geometry.Point(5.389, 52.077)),
('gb', shapely.geometry.Point(-1.131994, 52.866085)),
('us', shapely.geometry.Point(-102.743529, 42.077339)),
('an', shapely.geometry.Point(-68.87, 12.123)),
('ae', shapely.geometry.Point(54.469813, 23.549000 )),
('zw', shapely.geometry.Point(29.872, -19.0)),
('vn', shapely.geometry.Point(105.314, 21.491)),
('ua', shapely.geometry.Point(31.388, 49.016)),
('no', shapely.geometry.Point(8.74, 61.152)),
('li', shapely.geometry.Point(9.555, 47.153)),
('ma', shapely.geometry.Point(-5.758, 32.706)),
('es', shapely.geometry.Point(4.166025, 39.904935)),
('es', shapely.geometry.Point(-16.595450, 28.296571)),
('sd', shapely.geometry.Point(31.289617, 8.033712)),
('sd', shapely.geometry.Point(27.773992, 14.728712)),
('gb', shapely.geometry.Point(-6.811516, 58.000624)),
('gl', shapely.geometry.Point(-41.560707, 74.720801)),
('sm', shapely.geometry.Point(12.461467, 43.931073))
]
for code, point in points:
result = detect(point)
assert result['code'] == code
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment