diff -r 659ab465152e -r f9788970fa46 compute_regions.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/compute_regions.py	Wed Jul 29 23:45:53 2020 +0300
@@ -0,0 +1,99 @@
+#!/usr/bin/env python3
+import sys, json
+from misc import *
+from geometry import *
+from configparser import ConfigParser
+
+class Blockmap:
+	'''
+	The blockmap is a grid of latitude and longitude lines and models
+	a block -> set relation. A block is a latitude and longitude square.
+	'''
+	block_size = 200.0 # How big are blocks?
+	def __init__(self, blocks = None):
+		from collections import defaultdict
+		self.blocks = blocks or defaultdict(set)
+	def __getitem__(self, blockid):
+		'''
+		Returns a block for block coordinates. The block is a set that can
+		contain anything.
+		'''
+		return self.blocks[blockid]
+	def blockpoint(self, point):
+		'''
+		Returns blockmap coordinates for geographical coordinates.
+		The blockmap coordinates refer to a block in the blockmap.
+		'''
+		block_coordinate = lambda x: int(x * self.block_size)
+		return block_coordinate(point.x), block_coordinate(point.y)
+
+def blocks_in_shape(blockmap, shape):
+	'''
+	Finds all blocks inside the bounding box of a shape.
+	'''
+	from itertools import product
+	min_x, max_x = minmax(point.x for point in shape.points)
+	min_y, max_y = minmax(point.y for point in shape.points)
+	min_blockpoint = blockmap.blockpoint(Location(min_x, min_y))
+	max_blockpoint = blockmap.blockpoint(Location(max_x, max_y))
+	range_x = range(min_blockpoint[0], max_blockpoint[0] + 1)
+	range_y = range(min_blockpoint[1], max_blockpoint[1] + 1)
+	yield from (blockmap[x, y] for x, y in product(range_x, range_y))
+
+def create_blockmap(regions):
+	'''
+	Creates a blockmap of regions
+	'''
+	blockmap = Blockmap()
+	for region in regions.values():
+		# Minor shapes contain major shapes, so just use those
+		for shape in (region['minor_shapes'] or region['major_shapes']):
+			for block in blocks_in_shape(blockmap, shape):
+				set.add(block, region['name'])
+	return blockmap
+
+def get_args():
+	from argparse import ArgumentParser
+	parser = ArgumentParser()
+	parser.add_argument('gtfs_zip')
+	parser.add_argument('profile')
+	return parser.parse_args()
+
+def test_shapes(shapes, point):
+	return any(shape.contains_point(point) for shape in shapes)
+
+class RegionalLocation:
+	def __init__(self, *, region, region_class):
+		self.region, self.region_class = region, region_class
+	def __repr__(self):
+		return str.format(
+			'RegionalLocation(region = {region}, region_class = {region_class})',
+			region = repr(self.region),
+			region_class = repr(self.region_class),
+		)
+
+def locate_regionally(position, region):
+	if test_shapes(region['major_shapes'], position):
+		return RegionalLocation(region = region['name'], region_class = 'major')
+	elif test_shapes(region['minor_shapes'], position):
+		return RegionalLocation(region = region['name'], region_class = 'minor')
+	else:
+		return None
+
+def find_region_for_point(position, regions, blockmap):
+	for region_name in blockmap[blockmap.blockpoint(position)]:
+		region = regions[region_name]
+		classification = locate_regionally(position, region)
+		if classification:
+			return classification
+
+class RegionTester:
+	def __init__(self, regions):
+		self.regions = regions
+		self.blockmap = create_blockmap(regions)
+	def __call__(self, latitude, longitude):
+		return find_region_for_point(
+			position = Location(float(latitude), float(longitude)),
+			regions = self.regions,
+			blockmap = self.blockmap,
+		)