{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"from functools import partial\nimport time\nimport json\nimport os\nfrom glob import glob\nimport csv\nimport random\nimport urllib3\nimport itertools\n\nimport geopandas as gpd\nimport geojson\nimport folium\nimport requests\nimport requests_cache\nimport matplotlib as mpl\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport shapely\nfrom retry import retry\nfrom tqdm.autonotebook import tqdm\nimport pyproj\nimport shapely\nfrom functools import partial\nfrom shapely.geometry import box, Polygon, MultiPolygon, GeometryCollection, shape\nfrom shapely import wkt","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"cell_type":"code","source":"mpl.rcParams['figure.dpi'] = 200\npd.set_option('display.float_format', lambda x: f'{x:,.3f}')\npd.set_option(\"display.max_rows\", 100)\nurllib3.disable_warnings()\npd.set_option('display.width', 2000)\npd.set_option('max_colwidth', 400)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def point_to_latlon(point):\n    return point.y, point.x\n\ndef point_to_lonlat(point):\n    return point.x, point.y\n\ndef point_to_xy(point):\n    return point_to_lonlat(point)\n\ndef buffer_meters(geom, buffer):\n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={geom.centroid.y} +lon_0={geom.centroid.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    aeqd_to_wgs84 = partial(\n        pyproj.transform,\n        pyproj.Proj(local_azimuthal_projection),\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n    )\n\n    geom_transformed = shapely.ops.transform(wgs84_to_aeqd, geom)\n\n    buffer = geom_transformed.buffer(buffer)\n\n    buffer_wgs84 = shapely.ops.transform(aeqd_to_wgs84, buffer)\n    \n    return buffer_wgs84  ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_city = gpd.read_file('../data/city_boundaries_100m.geojson')\ndf_city.sort_values('area_km2').tail(10)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"iloc = 0\npoly = df_city.iloc[iloc].geometry\n\n\n# Plot point and geom.\nfig = folium.Figure(width=1000, height=500)\nfm = folium.Map(\n    tiles='OpenStreetMap',\n    location=point_to_latlon(poly.centroid),\n    zoom_start=10,\n).add_to(fig)\n\nfolium.GeoJson(poly).add_to(fm)\n\n\nfm","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Adapted from https://github.com/JoaoCarabetta/osm-road-length\n\ndef simplify(geom):\n    local_azimuthal_projection = f\"+proj=aeqd +R=6371000 +units=m +lat_0={geom.centroid.y} +lon_0={geom.centroid.x}\"\n\n    wgs84_to_aeqd = partial(\n        pyproj.transform,\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n        pyproj.Proj(local_azimuthal_projection),\n    )\n\n    aeqd_to_wgs84 = partial(\n        pyproj.transform,\n        pyproj.Proj(local_azimuthal_projection),\n        pyproj.Proj('+proj=longlat +datum=WGS84 +no_defs'),\n    )\n\n    geom_transformed = shapely.ops.transform(wgs84_to_aeqd, geom)\n    geom_simple = geom_transformed.simplify(500)\n    geom_wgs84 = shapely.ops.transform(aeqd_to_wgs84, geom_simple)\n    \n    return geom_wgs84  \n\n\n\n\ndef to_geojson(x):\n\n    if isinstance(x, shapely.geometry.multipolygon.MultiPolygon):\n        x = max(x, key=lambda a: a.area)\n\n    g = geojson.Feature(geometry=x, properties={}).geometry\n\n    return g[\"coordinates\"][0]\n\n\ndef swipe(x):\n    return [[c[1], c[0]] for c in x]\n\n\ndef flatten(l):\n    return [str(round(item, 4)) for sublist in l for item in sublist]\n\n\ndef to_overpass_coords(x):\n\n    coords = to_geojson(x)\n    coords = swipe(coords)\n    coords = flatten(coords)\n    coords = \" \".join(coords)\n    return coords\n\ndef threshold_func(g, value):\n\n    return get_area(g) < value\n\ndef katana(geometry, threshold_func, threshold_value, count=0):\n    \"\"\"Split a Polygon into two parts across it's shortest dimension\n    \n    KUDOS https://snorfalorpagus.net/blog/2016/03/13/splitting-large-polygons-for-faster-intersections/\n    \"\"\"\n    bounds = geometry.bounds\n    width = bounds[2] - bounds[0]\n    height = bounds[3] - bounds[1]\n    if threshold_func(geometry, threshold_value) or count == 250:\n        # either the polygon is smaller than the threshold, or the maximum\n        # number of recursions has been reached\n        return [geometry]\n    if height >= width:\n        # split left to right\n        a = box(bounds[0], bounds[1], bounds[2], bounds[1] + height / 2)\n        b = box(bounds[0], bounds[1] + height / 2, bounds[2], bounds[3])\n    else:\n        # split top to bottom\n        a = box(bounds[0], bounds[1], bounds[0] + width / 2, bounds[3])\n        b = box(bounds[0] + width / 2, bounds[1], bounds[2], bounds[3])\n    result = []\n    for d in (\n        a,\n        b,\n    ):\n        c = geometry.intersection(d)\n        if not isinstance(c, GeometryCollection):\n            c = [c]\n        for e in c:\n            if isinstance(e, (Polygon, MultiPolygon)):\n                result.extend(katana(e, threshold_func, threshold_value, count + 1))\n    if count > 0:\n        return result\n    # convert multipart into singlepart\n    final_result = []\n    for g in result:\n        if isinstance(g, MultiPolygon):\n            final_result.extend(g)\n        else:\n            final_result.append(g)\n    return final_result\n\n\n\ndef get_area(s):\n    # s = shape(s)\n    proj = partial(\n        pyproj.transform, pyproj.Proj(\"epsg:4326\"), pyproj.Proj(\"epsg:3857\")\n    )\n\n    return shapely.ops.transform(proj, s).area / 1e6  # km\n\n\n@retry(tries=2, delay=1, backoff=2, max_delay=60)\ndef overpass_request(geo, yomo=None, overpass_url='http://overpass-api.de/api/interpreter'):\n    \n    if yomo:\n        date_str = f'[date:\"{yomo}-01T00:00:00Z\"]'\n    else:\n        date_str = ''\n        \n    overpass_query = f\"\"\"\n        [out:json][timeout:25];\n        (\n          node[\"amenity\"=\"hospital\"][\"emergency\"=\"yes\"](poly:\"{geo}\");\n          way[\"amenity\"=\"hospital\"][\"emergency\"=\"yes\"](poly:\"{geo}\");\n          relation[\"amenity\"=\"hospital\"][\"emercency\"=\"yes\"](poly:\"{geo}\");\n        );\n        out body;\n        >;\n        out skel qt;\n    \"\"\"\n        \n\n   \n    \n    response = requests.get(\n        overpass_url,\n        params={\"data\": overpass_query},\n        proxies=proxies,\n        verify=False,\n        timeout=60,\n    )\n    \n    return response\n\n\ndef get(geometry, yomo=None, overpass_url='http://overpass-api.de/api/interpreter', threshold_value=1000000/100):\n    \"\"\"Get Open Street Maps highways length in meters and object count for a given geometry\n    It splits the regions to manage overpass turbo limits.\n    For MultiPolygons, only the biggest polygon will be considered.\n    Parameters\n    ----------\n    geometry : shapely.geometry\n        A shapely polygon\n    threshold_value : int, optional\n        Maximum area in sq km to split the polygons, by default 1000000\n    Returns\n    -------\n    pd.DataFrame\n        Table indexed by highway with length sum in meters and observation count\n    \"\"\"\n\n    geo = simplify(geometry)\n    geoms = katana(geometry, threshold_func, threshold_value)\n    responses = []\n    print(len(geoms))\n    for geo in geoms:\n        geo = to_overpass_coords(geo)\n        response = overpass_request(geo, yomo, overpass_url)\n        responses.append([])\n        \n    \n    return responses\n    \n    \n# r = get(poly)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ilocs = list(range(len(df_city)))\nrandom.shuffle(ilocs)\n        \nbackends = [\n    'http://lz4.overpass-api.de/api/interpreter', \n    'http://z.overpass-api.de/api/interpreter',\n    'https://overpass.openstreetmap.ru/api/interpreter',\n    'https://overpass.nchc.org.tw/api/interpreter',\n]\nbackends_r = itertools.cycle(backends)\n        \n    \nfor iloc in tqdm(ilocs):\n    path = f'../data/hosp/{iloc}.csv'\n    if os.path.exists(path):\n        continue\n    overpass_url = next(backends_r)\n    poly = buffer_meters(df_city.iloc[iloc].geometry, 500)\n\n    responses = get(poly, overpass_url=overpass_url)\n\n    dfs = [pd.DataFrame.from_records(response.json()['elements']) for response in responses]\n    df = pd.concat(dfs)\n    df['city_iloc'] = iloc\n\n    df.to_csv(path, index=False, quoting=csv.QUOTE_NONNUMERIC)\n    time.sleep(10)\n        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Merge into one df.\ndfs = []\nfor iloc in tqdm(ilocs):\n    path = f'../data/hosp/{iloc}.csv'\n    if not os.path.exists(path):\n        continue\n        \n    df = pd.read_csv(path)\n    if len(df) == 0:\n        continue\n    dfs.append(df)\n    \ndf_merged = pd.concat(dfs)\ndf_merged.to_csv('../data/osrm_hospital_location.csv')","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}