Nous publierons un programme qui recherche des polygones contenant certaines coordonnées en utilisant les données de polygones publiées. Obtenir une adresse à partir de la latitude et de la longitude s'appelle le géocodage inversé.
Si vous ne disposez pas de données polygonales, le géocodage inversé est également fourni par l'API Google Maps, etc. et il est facile à utiliser. Si vous disposez de données polygonales et que vous souhaitez effectuer un géocodage inversé, vous pouvez également utiliser la fonction SIG d'une base de données telle que MySQL.
Cette fois, il s'agit d'un programme Python lorsque vous souhaitez déterminer quelle ville, quartier, ville ou village au Japon appartient à un grand nombre de coordonnées. Même un petit programme fish qui génère un index à chaque fois a beaucoup de coordonnées, donc il a fallu moins de temps pour s'exécuter que lors de l'utilisation de MySQL 5.6.
Avantages de l'utilisation de MySQL:
Démérite:
――Vous avez besoin de données de polygones qui correspondent au niveau d'adresse souhaité (les polygones de niveau d'adresse peuvent ne pas être facilement disponibles)
Raison de son implémentation en Python:
――Vous n'avez pas besoin d'être aussi détaillé que le niveau d'adresse --Lorsque vous accédez à la base de données, le montant est trop important pour terminer (bien sûr, l'API est également difficile) --La source des données de zone utilisées pour le classement est claire
Pip install rtree shapely` `` seul peut ne pas fonctionner car il y a des bibliothèques dépendantes. Veuillez donc lire chaque document pour l'installation.
Vous pouvez écrire ceci en utilisant rtree et galbé.
revgeocoder.py
# -*- coding: utf-8 -*-
import collections
from shapely.geometry import Polygon, Point
from rtree import index
Area = collections.namedtuple('Area', ['area_id', 'polygon'])
class ReverseGeocoder():
def __init__(self):
self.idx = index.Index()
def insert_from_iterator(self, itr):
'''(id, Polygon)Créer un Rtree à partir d'un itérateur qui retourne
Polygon.Rtree créé en fonction de l'id et du polygone liés.
Args:
itr: (id, Polygon)Itérateur qui retourne
'''
for i, (area_id, polygon) in enumerate(itr):
obj = Area(area_id=area_id, polygon=polygon)
self.idx.insert(i, polygon.bounds, obj)
def contains(self, lat, lon):
'''Point(lat, lon)Zone de polygone comprenant_Renvoie id.
S'il correspond à deux polygones ou plus, il est trié et renvoyé par ordre croissant de superficie.
(Lorsqu'une zone a une relation d'inclusion, il est préférable de sélectionner celle avec la zone la plus petite.
Cependant, en réalité, il y a un chevauchement entre Polygon et Polygon ... )
'''
result = []
point = Point(lat, lon)
for hit in self.idx.intersection(point.bounds, objects=True):
if hit.object.polygon.contains(point):
result.append(hit.object)
if len(result) > 1:
result.sort(key=lambda x: (x.polygon.area, x.area_id))
return [r.area_id for r in result]
def __repr__(self):
return '<ReverseGeocoder contains {} polygons>'.format(self.idx.count(self.idx.bounds))
Lorsque vous utilisez les données de limite qui peuvent être téléchargées depuis e-stat, cela ressemble à ceci.
findcity.py
# -*- coding: utf-8 -*-
import glob
import fiona
from shapely.geometry import Polygon, shape
from revgeocoder import ReverseGeocoder
def parse_shapefile(shapefile):
'''(area_id, Polygon, name)Générateur qui retourne'''
with fiona.open(shapefile, 'r') as source:
for obj in source:
#Données polygonales
polygon = shape(obj['geometry'])
#Adresse de la zone
names = []
for prop in ['KEN_NAME', 'GST_NAME', 'CSS_NAME', 'MOJI']:
if obj['properties'][prop] is not None:
names.append(obj['properties'][prop])
name = ''.join(names)
#Générer un identifiant unique (KEN+ CITY + SEQ_NO2）
area_id = int(''.join(map(str, (obj['properties']['KEN'], obj['properties']['CITY'], obj['properties']['SEQ_NO2']))))
yield area_id, polygon, name
def main():
shapefiles = glob.glob('data/japan-shape/A002005212010DDSWC3520*/*.shp')
print('Shapefiles:', shapefiles)
#Faire un index
rgeocoder = ReverseGeocoder()
id2name = {}
def gen(shapefile):
for area_id, polygon, name in parse_shapefile(shapefile):
id2name[area_id] = name
yield area_id, polygon
for shapefile in shapefiles:
rgeocoder.insert_from_iterator(gen(shapefile))
# test
area_id = rgeocoder.contains(132.257269, 34.108815)[0]
print(area_id, id2name[area_id])
return rgeocoder
if __name__ == '__main__':
rgeocoder = main()
Exemple d'exécution:
$ python findcity.py
Shapefiles: ['data/japan-shape/A002005212010DDSWC35208/h22ka35208.shp', 'data/japan-shape/A002005212010DDSWC35207/h22ka35207.shp', 'data/japan-shape/A002005212010DDSWC35201/h22ka35201.shp', 'data/japan-shape/A002005212010DDSWC35203/h22ka35203.shp', 'data/japan-shape/A002005212010DDSWC35206/h22ka35206.shp', 'data/japan-shape/A002005212010DDSWC35202/h22ka35202.shp', 'data/japan-shape/A002005212010DDSWC35204/h22ka35204.shp']
35208602 6-chome, Sofucho, ville d'Iwakuni, préfecture de Yamaguchi
Recommended Posts