{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import copy\nimport json\nimport re\nimport time\nimport unicodedata\nimport warnings\nfrom functools import partial\nimport os\nimport math\n\nfrom tqdm.autonotebook import tqdm\nfrom urllib.parse import urlencode\nimport contextily\nimport folium\nimport geopandas as gpd\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport pyproj\nimport rasterio\nimport rasterio.warp\nimport rasterio.features\nimport rasterio.transform\nimport rasterio.plot\nimport requests\nimport requests_cache\nimport seaborn as sns\nimport shapely.wkt\nimport polyline\nimport scipy.interpolate\n\n%matplotlib inline","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":false},"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 \n\n\ndef lat_lon_distance(from_lat, from_lon, to_lat, to_lon):\n    '''\n    Find the distance between two latlon points.\n    Uses the Haversine forumula.\n    to_lat and to_lon can be vectors of points.\n    '''\n    mean_lat = (from_lat + to_lat) / 2\n    radius = radius_of_earth(mean_lat)\n\n    from_lat = math.radians(from_lat)\n    from_lon = math.radians(from_lon)\n    to_lat = np.radians(to_lat)\n    to_lon = np.radians(to_lon)\n    d_lat = to_lat - from_lat\n    d_lon = to_lon - from_lon\n\n    a = np.sin(d_lat / 2)**2 + math.cos(from_lat) * np.cos(to_lat) * np.sin(d_lon  / 2)**2\n    d = radius * 2 * np.arcsin(np.sqrt(a))\n    return d\n\n\ndef radius_of_earth(lat):\n    '''\n    Find the radius of the Earth at a given latitude. Vectorised.\n    Formula and data from https://en.wikipedia.org/wiki/Earth_radius\n    '''\n    r_1 = 6378137.0  # Earth equatorial radius\n    r_2 = 6356752.3  # Earth polar radius\n\n    lat = np.radians(lat)\n    part_1 = r_1**4 * np.cos(lat)**2 + r_2**4 * np.sin(lat)**2\n    part_2 = r_1**2 * np.cos(lat)**2 + r_2**2 * np.sin(lat)**2\n\n    return np.sqrt(part_1 / part_2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_city = gpd.read_file('../data/city_boundaries_100m.geojson')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_hospitals(iloc):\n    path = f'../data/hosp/{iloc}.csv'\n    if not os.path.exists(path):\n        return None\n    \n    with open(path) as f:\n        df_hosp = pd.read_csv(path)\n    \n    \n    hospitals = []\n    for t in df_hosp.itertuples():\n        if t.type == 'node' and not pd.isnull(t.lat) and not pd.isnull(t.lon) and not pd.isnull(t.tags):\n            hospitals.append((t.lat, t.lon))\n            continue\n\n        if t.type == 'way':\n            node_id = eval(t.nodes)[0]\n            node = df_hosp[df_hosp.id == node_id].iloc[0]\n            hospitals.append((node.lat, node.lon))\n            \n    # Remove too close together.\n    lats = np.array([h[0] for h in hospitals])\n    lons = np.array([h[1] for h in hospitals])\n    cluster_threshold = 500\n\n    remove_point = np.full(len(lats), False, dtype=bool)\n    for i, (lat, lon) in enumerate(zip(lats, lons)):\n        if remove_point[i]:\n            continue\n        distances = lat_lon_distance(lat, lon, lats[i+1:], lons[i+1:])\n        remove_point[i+1:] = remove_point[i+1:] | (distances < cluster_threshold)\n\n    lats = lats[~remove_point]\n    lons = lons[~remove_point]\n    hospitals = list(zip(lats, lons))\n    \n            \n    return hospitals\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"POP_PATH_1k = '../input/ppp_2020_1km_Aggregated.tif'\nPOP_PATH_100m = '../input/MOSAIC_ppp_prj_2020/MOSAIC_ppp_prj_2020.vrt'\ndef load_pop_raster(geom, buffer=1_000, pop_path=POP_PATH_100m):\n\n    # Load population raster.\n    geom_buffer = buffer_meters(geom, buffer)\n    with rasterio.open(pop_path) as f:\n        window = rasterio.windows.from_bounds(*geom_buffer.bounds, f.transform)\n        a = f.read(1, window=window, masked=True)\n        a = np.ma.filled(a, 0)\n        \n\n    return f, a, window\n\ndef load_mortality_raster(geom, buffer=1_000):\n    path = '../input/povmap-global-subnational-infant-mortality-rates-v2-geotiff/povmap_global_subnational_infant_mortality_rates_v2.tif'\n    geom_buffer = buffer_meters(geom, buffer)\n    with rasterio.open(path) as f:\n        window = rasterio.windows.from_bounds(*geom_buffer.bounds, f.transform)\n        a = f.read(1, window=window, masked=True)\n        a = np.ma.filled(a, 0)\n        a[a < 0] = 0\n        \n    return f, a, window\n\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        \n        \n       \ndef weighted_raster_sample(f, a_geom, window, n_samples):\n    weights = a_geom / a_geom.sum()\n    weights_flat = weights.flatten()\n    lonlats = []\n    n = len(weights_flat)\n    for _ in range(n_samples):\n        idx_flat = np.random.choice(np.arange(n), p=weights_flat)\n        idx = np.unravel_index(idx_flat, a_geom.shape)\n\n        # Add subcell jitter.\n        j0 = np.random.uniform(0, 1)\n        j1 = np.random.uniform(0, 1)\n        rowcol = (idx[0] + j0, idx[1] + j1)\n        lonlat = rasterio.transform.xy(f.window_transform(window), rowcol[0], rowcol[1])\n        lonlats.append(lonlat)\n        \n    return lonlats\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Pop-weighted sample.\n\nfor iloc in tqdm(range(len(df_city))):\n    osrm_path = f'../data/osrm-low/{iloc}.csv'\n    if os.path.exists(osrm_path):\n        continue\n        \n    # Load hospitals.\n    hospitals = load_hospitals(iloc)\n    if len(hospitals) == 0:\n        continue\n    \n    # Load population density raster.\n    geom = df_city.geometry.iloc[iloc]\n    f, a, window = load_pop_raster(geom)\n    if max(a.shape) < 100:\n        f, a, window = load_pop_raster(geom, pop_path=POP_PATH_100m)\n    mask = rasterio.features.geometry_mask([geom], a.shape, f.window_transform(window), invert=True)\n    a_geom = a * mask\n    \n    # Load infant mortality raster.\n    f_im, a_im, window_im = load_mortality_raster(geom)\n    x_im = np.arange(a_im.shape[1])\n    y_im = np.arange(a_im.shape[0])\n    x_a = np.arange(a.shape[1])\n    y_a = np.arange(a.shape[0])\n    x_im = np.linspace(0, 1, a_im.shape[1])\n    y_im = np.linspace(0, 1, a_im.shape[0])\n    x_a = np.linspace(0, 1, a.shape[1])\n    y_a = np.linspace(0, 1, a.shape[0])\n    a_im_intep = scipy.interpolate.RectBivariateSpline(y_im, x_im, a_im, kx=1, ky=1)(y_a, x_a)\n    a_im_intep = a_im_intep + np.random.uniform(0, a_im_intep.max()/100, a_im_intep.shape)\n    a_low_im = (a_im_intep < np.median(a_im_intep[mask])).astype(int)\n    a_geom = a_geom * a_low_im\n    \n    # Pop-weighted random sample.\n    source_lonlats = weighted_raster_sample(f, a_geom, window, n_samples=250)\n    \n    # Build OSM query.\n    source_lats = [p[1] for p in source_lonlats]\n    source_lons = [p[0] for p in source_lonlats]\n    dest_lats = [p[0] for p in hospitals]\n    dest_lons = [p[1] for p in hospitals]\n    latlons = [(lat, lon) for lon, lat in source_lonlats] + hospitals\n    pline = polyline.encode(latlons)\n    sources = ';'.join(str(i) for i in range(len(source_lats)))\n    destinations = ';'.join(str(i + len(source_lats)) for i in range(len(dest_lats)))\n    url = f'http://router.project-osrm.org/table/v1/driving/polyline({pline})?sources={sources}&destinations={destinations}'\n    \n    response = requests.get(url)\n    time.sleep(10)\n    \n    # Save dataframe.\n    records = []\n    durations = np.array(response.json()['durations'])\n    for i_source in range(len(source_lats)):\n        for i_dest in range(len(dest_lats)):\n            records.append({\n                'source_lat': source_lats[i_source],\n                'source_lon': source_lons[i_source],\n                'dest_lat': dest_lats[i_dest],\n                'dest_lon': dest_lons[i_dest],\n                'duration': durations[i_source, i_dest]\n            })\n    df_osrm = pd.DataFrame.from_records(records)\n    df_osrm.to_csv(osrm_path, index=False)\n        \n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Merge into one table.\ndfs = []\nfor iloc in tqdm(range(len(df_city))):\n    osrm_path = f'../data/osrm-low/{iloc}.csv'\n    if not os.path.exists(osrm_path):\n        continue\n        \n    df = pd.read_csv(osrm_path)\n    df['city_iloc'] = iloc\n    dfs.append(df)\n        \ndf_merged = pd.concat(dfs)\ndf_merged.to_csv('../data/osrm_hospital_duration_low_im.csv', index=False)","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}