We will publish a program that searches for polygons containing certain coordinates using the published polygon data. Getting an address from latitude and longitude is called Reverse geocoding.
If you don't have polygon data, Reverse geocoding is also provided by Google Maps API etc. and it is easy to use. If you have polygon data and want to reverse geocode, you can also use the GIS function of a database such as MySQL.
This time, it is a Python program when you want to determine which city, ward, town, or village in Japan belongs to a large number of coordinates. Even a small fish program that generates an index every time has many coordinates, so it took less time to execute than when using MySQL 5.6.
Benefits of using MySQL:
--If you have data, just put the data and write SQL
Demerit:
--You need polygon data that matches the address level you want (address level polygons may not be available easily)
Why implemented in Python:
――You don't have to be as detailed as the address level --When you access the database, the amount is too large to finish (of course, API is also difficult) --The source of the area data used for verification is clear
`pip install rtree shapely` alone may not work because there are dependent libraries. So please read each document for installation.
You can write this using rtree and shapely.
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)Create a Rtree from an iterator that returns
Polygon.Rtree created based on bound holds id and Polygon.
Args:
itr: (id, Polygon)Iterator that returns
'''
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)Polygon area including_Returns id.
If it matches two or more Polygons, it is sorted and returned in ascending order of area.
(When there is an area that has an inclusion relationship, it is better to select the one with the smaller area.
However, in reality, there is an overlap between Polygon and 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))
This is the case when using boundary data that can be downloaded from e-stat.
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)Generator that returns'''
with fiona.open(shapefile, 'r') as source:
for obj in source:
#Polygon data
polygon = shape(obj['geometry'])
#Address of the area
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)
#Generate a unique id (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)
#Make an 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()
Execution example:
$ 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, Shozokumachi, Iwakuni City, Yamaguchi Prefecture