Created
March 16, 2014 11:28
-
-
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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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