terrassumptions/server/app/point_gen/urban_centers.py

import math
import random
import csv
import logging
from collections import defaultdict

from .shared import point_has_streetview
from ..scoring import mean_earth_radius_km

logger = logging.getLogger(__name__)

URBAN_CENTERS = defaultdict(list)
_found_countries = set()
_urban_center_count = 0
with open("./data/urban-centers.csv") as infile:
    for code, name, lat, lng in csv.reader(infile, delimiter=",", quotechar='"'):
        URBAN_CENTERS[code].append((name, float(lat), float(lng)))
        _found_countries.add(code)
        _urban_center_count += 1
logger.info(f"Read {_urban_center_count} urban centers from {len(_found_countries)} countries.")
VALID_COUNTRIES = tuple(_found_countries)


async def urban_coord(country_lock, city_retries=10, point_retries=10, max_dist_km=25):
    """
    Returns (latitude, longitude) of usable coord (where google has data) that is near
    a known urban center. Points will be at most max_dist_km kilometers away. This function 
    will use country_lock to determine the country from which to pull a known urban center,
    generate at most point_retries points around that urban center, and try at most
    city_retries urban centers in that country. If none of the generated points have street 
    view data, this will return None. Otherwise, it will exit as soon as suitable point is 
    found.

    This function calls the streetview metadata endpoint - there is no quota consumed.
    """

    country_lock = country_lock.lower()
    cities = URBAN_CENTERS[country_lock]
    src = random.sample(cities, k=min(city_retries, len(cities)))

    logger.info(f"Trying {len(src)} centers in {country_lock}")

    for (name, city_lat, city_lng) in src:
        # logic adapted from https://stackoverflow.com/a/7835325
        # start in a city
        logger.info(f"Trying at most {point_retries} points around {name}")
        city_lat_rad = math.radians(city_lat)
        sin_lat = math.sin(city_lat_rad)
        cos_lat = math.cos(city_lat_rad)
        city_lng_rad = math.radians(city_lng)
        for _ in range(point_retries):
            # turn a random direction, and go random distance
            dist_km = random.random() * max_dist_km
            angle_rad = random.random() * 2 * math.pi
            d_over_radius = dist_km / mean_earth_radius_km
            sin_dor = math.sin(d_over_radius)
            cos_dor = math.cos(d_over_radius)
            pt_lat_rad = math.asin(sin_lat * cos_dor + cos_lat * sin_dor * math.cos(angle_rad))
            pt_lng_rad = city_lng_rad + math.atan2(math.sin(angle_rad) * sin_dor * cos_lat, cos_dor - sin_lat * math.sin(pt_lat_rad))
            pt_lat = math.degrees(pt_lat_rad)
            pt_lng = math.degrees(pt_lng_rad)
            if await point_has_streetview(pt_lat, pt_lng):
                logger.info("Point found!")
                return (country_lock, pt_lat, pt_lng)