diff --git a/notebooks/IMPLEMENTATIONS/TUN_Road_Improvements/Friction_Generation.ipynb b/notebooks/IMPLEMENTATIONS/TUN_Road_Improvements/Friction_Generation.ipynb new file mode 100644 index 0000000..80fa591 --- /dev/null +++ b/notebooks/IMPLEMENTATIONS/TUN_Road_Improvements/Friction_Generation.ipynb @@ -0,0 +1,442 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculating travel time to ports based on select road improvements in Tunisia\n", + "\n", + "Unit of analysis is ADM2 (geoboundaries). Need to summarize the following:\n", + "- Monthly Nighttime Lights mean and sum\n", + "- TT to ports pre and post implementation of TREATED roads\n", + "- Distance to nearest treated road, distance to each labelled treated road\n", + "- Percentage Urban" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import sys, os, importlib\n", + "import rasterio\n", + "\n", + "import pandas as pd\n", + "import geopandas as gpd\n", + "import skimage.graph as graph\n", + "import numpy as np\n", + "\n", + "from shapely.geometry import Point\n", + "\n", + "sys.path.insert(0, r\"C:\\WBG\\Work\\Code\\GOSTrocks\\src\")\n", + "import GOSTrocks.rasterMisc as rMisc\n", + "import GOSTrocks.osmMisc as osmMisc\n", + "import GOSTrocks.dataMisc as dMisc\n", + "import GOSTrocks.ntlMisc as ntlMisc\n", + "from GOSTrocks.misc import tPrint\n", + "\n", + "sys.path.append(r\"C:\\WBG\\Work\\Code\\GOSTnetsraster\\src\")\n", + "import GOSTnetsraster.market_access as ma\n", + "import GOSTnetsraster.conversion_tables as speed_tables\n", + "\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Input parameters\n", + "m_crs = 32632 # Need to project data to a metres-based projection\n", + "\n", + "# Define input data\n", + "base_folder = \"C:/WBG/Work/TUN_Impact/\"\n", + "landcover_file = os.path.join(base_folder, \"DATA\", 'ESA_Globcover.tif')\n", + "# These are the digitized road segements that have been improved\n", + "treated_segments_file = os.path.join(base_folder, \"DATA\", \"DIME_Roads\", 'treated_roads.shp')\n", + "control_segments_file = os.path.join(base_folder, \"DATA\", \"DIME_Roads\", 'control_roads.shp')\n", + "road_segments_file = os.path.join(base_folder, \"DATA\", \"impacted_osm_roads.gpkg\")\n", + "osm_roads_file = os.path.join(base_folder, \"DATA\", \"OSM\", \"gis_osm_roads_free_1.shp\")\n", + "# WorldPop 2020 constrained, projected to m_crs\n", + "pop_file = os.path.join(base_folder, \"DATA\", \"tun_ppp_2020_constrained_proj.tif\")\n", + "# https://datacatalog.worldbank.org/int/search/dataset/0038118/Global---International-Ports\n", + "port_file = os.path.join(base_folder, \"DATA\", \"TUN_ports.gpkg\")\n", + "# administrative bounadaries are used to summarize population\n", + "tun_adm2 = dMisc.get_geoboundaries(\"TUN\", 'ADM2')\n", + "tun_adm1 = dMisc.get_geoboundaries(\"TUN\", 'ADM1')\n", + "\n", + "# Define output files\n", + "friction_folder = os.path.join(base_folder, \"DATA\", \"FRICTION\")\n", + "results_folder = os.path.join(base_folder, \"RESULTS\")\n", + "for cFolder in [friction_folder, results_folder]:\n", + " if not os.path.exists(cFolder):\n", + " os.makedirs(cFolder) \n", + "pre_friction_file = os.path.join(friction_folder, 'FRICTION_pre_intervention.tif')\n", + "post_friction_file = os.path.join(friction_folder, 'FRICTION_post_intervention.tif')\n", + "# This extracts the existing global friction file, used only for comparison\n", + "global_friction_file = os.path.join(friction_folder, \"2020_motorized_friction.geotiff\")\n", + "if not os.path.exists(global_friction_file):\n", + " gl_fr = rasterio.open(r\"J:\\Data\\GLOBAL\\INFRA\\FRICTION_2020\\2020_motorized_friction_surface.geotiff\")\n", + " local_fr = rMisc.clipRaster(gl_fr, tun_adm2, global_friction_file)\n", + " \n", + "# Read in data\n", + "dests = gpd.read_file(port_file).to_crs(m_crs)\n", + "if not os.path.exists(landcover_file):\n", + " global_landcover = r\"R:\\GLOBAL\\LCVR\\Globcover\\2015\\ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif\"\n", + " in_lc = rasterio.open(global_landcover)\n", + " temp_landcover_file = landcover_file.replace(\".tif\", \"_temp.tif\")\n", + " local_lc = rMisc.clipRaster(in_lc, tun_adm2, temp_landcover_file)\n", + " temp_lc = rasterio.open(temp_landcover_file)\n", + " proj_res = rMisc.project_raster(temp_lc, m_crs)\n", + " with rasterio.open(landcover_file, 'w', **proj_res[1]) as outR:\n", + " outR.write(proj_res[0])\n", + "\n", + "in_lc = rasterio.open(landcover_file)\n", + "in_pop = rasterio.open(pop_file)\n", + "if in_pop.crs != in_lc.crs:\n", + " proj_res = rMisc.standardizeInputRasters(in_pop, in_lc, pop_file.replace(\".tif\", \"_proj.tif\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Downlaod worldcover data\n", + "tiles_geojson = r\"C:\\WBG\\Work\\data\\LCVR\\esa_worldcover_grid.geojson\"\n", + "in_tiles = gpd.read_file(tiles_geojson)\n", + "sel_tiles = in_tiles.loc[in_tiles.intersects(tun_adm2.unary_union)]\n", + "\n", + "tile_path = \"s3://esa-worldcover/v200/2021/map/ESA_WorldCover_10m_2021_v200_{tile}_Map.tif\"\n", + "out_folder = os.path.join(base_folder, \"DATA\", \"WorldCover\")\n", + "for idx, row in sel_tiles.iterrows():\n", + " cur_tile_path = tile_path.format(tile=row['ll_tile'])\n", + " cur_out = os.path.join(out_folder, f\"WorldCover_{row['ll_tile']}.tif\")\n", + " if not os.path.exists(cur_out):\n", + " command = f\"aws s3 --no-sign-request --no-verify-ssl cp {cur_tile_path} {cur_out}\"\n", + " print(command)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Download and process industrial zones\n", + "import urllib.request, json \n", + "with urllib.request.urlopen(\"https://afi.e-industrie.gov.tn/apps-lots.php\") as url:\n", + " industry_zones = json.load(url)\n", + "zones_df = pd.DataFrame(industry_zones)\n", + "\n", + "zones_geoms = [Point(x) for x in zip(zones_df['ZONE_LONG'], zones_df['ZONE_LAT'])]\n", + "zones_df = gpd.GeoDataFrame(zones_df, geometry=zones_geoms, crs=4326)\n", + "zones_df = zones_df.to_crs(m_crs)\n", + "zones_df.to_file(os.path.join(base_folder, \"DATA\", \"MAPPING\", \"industrial_zones.gpkg\"), driver='GPKG')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Process roads to create pre and post friction surfaces\n", + "sel_roads = gpd.read_file(road_segments_file)\n", + "sel_roads = sel_roads.to_crs(m_crs)\n", + "sel_roads['speed'] = 10\n", + "all_roads = gpd.read_file(osm_roads_file)\n", + "all_roads = all_roads.to_crs(m_crs)\n", + "all_roads['speed'] = all_roads['fclass'].map(speed_tables.osm_speed_dict)\n", + "all_roads['speed'] = all_roads['speed'].fillna(10.0)\n", + "wb_roads_ids = sel_roads.loc[~sel_roads['osm_id'].isna(),'osm_id']\n", + "new_roads = sel_roads.loc[sel_roads['osm_id'].isna(),]\n", + "\n", + "lc_speed_table = speed_tables.esaacci_landcover" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Generate pre-intervention friction surface\n", + "if not os.path.exists(pre_friction_file):\n", + " pre_roads = all_roads.copy()\n", + " pre_roads.loc[pre_roads['osm_id'].isin(wb_roads_ids), 'speed'] = 10.0\n", + "\n", + " pre_friction = ma.generate_roads_lc_friction(in_lc, pre_roads, lc_travel_table=lc_speed_table, \n", + " out_file=pre_friction_file, resolution=in_lc.res[0])\n", + "\n", + "pre_friction = rasterio.open(pre_friction_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Generate post-intervention friction surface\n", + "if not os.path.exists(post_friction_file):\n", + " post_roads = all_roads.copy()\n", + " post_roads.loc[post_roads['osm_id'].isin(wb_roads_ids), 'speed'] = 40.0\n", + "\n", + " post_friction = ma.generate_roads_lc_friction(in_lc, post_roads, lc_travel_table=lc_speed_table, \n", + " out_file=post_friction_file, resolution=in_lc.res[0])\n", + " \n", + "post_friction = rasterio.open(post_friction_file)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculate travel time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate pre-intervention, population-weighted travel time\n", + "frictionD = pre_friction.read()[0,:,:]\n", + "frictionD = frictionD * pre_friction.res[0]\n", + "mcp = graph.MCP_Geometric(frictionD)\n", + "pre_tt_ports = ma.summarize_travel_time_populations(in_pop, pre_friction, dests, mcp, tun_adm2)\n", + "pd.DataFrame(pre_tt_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"PRE_ADM2_tt_ports.csv\"))\n", + "\n", + "pre_zones_ports = ma.summarize_travel_time_populations(in_pop, pre_friction, zones_df, mcp, tun_adm2)\n", + "pd.DataFrame(pre_zones_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"PRE_ADM2_tt_zones.csv\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate pre-intervention, population-weighted travel time\n", + "frictionD = pre_friction.read()[0,:,:]\n", + "frictionD = frictionD * pre_friction.res[0]\n", + "mcp = graph.MCP_Geometric(frictionD)\n", + "pre_tt_ports = ma.summarize_travel_time_populations(in_pop, pre_friction, dests, mcp, tun_adm1)\n", + "pd.DataFrame(pre_tt_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"PRE_ADM1_tt_ports.csv\"))\n", + "\n", + "pre_zones_ports = ma.summarize_travel_time_populations(in_pop, pre_friction, zones_df, mcp, tun_adm1)\n", + "pd.DataFrame(pre_zones_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"PRE_ADM1_tt_zones.csv\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate pre-intervention, population-weighted travel time\n", + "frictionD = post_friction.read()[0,:,:]\n", + "frictionD = frictionD * post_friction.res[0]\n", + "mcp = graph.MCP_Geometric(frictionD)\n", + "post_tt_ports = ma.summarize_travel_time_populations(in_pop, post_friction, dests, mcp, tun_adm2)\n", + "pd.DataFrame(post_tt_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"POST_ADM2_tt_ports.csv\"))\n", + "\n", + "post_zones_ports = ma.summarize_travel_time_populations(in_pop, post_friction, zones_df, mcp, tun_adm2)\n", + "pd.DataFrame(post_zones_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"POST_ADM2_tt_zones.csv\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate pre-intervention, population-weighted travel time\n", + "frictionD = post_friction.read()[0,:,:]\n", + "frictionD = frictionD * post_friction.res[0]\n", + "mcp = graph.MCP_Geometric(frictionD)\n", + "post_tt_ports = ma.summarize_travel_time_populations(in_pop, post_friction, dests, mcp, tun_adm1)\n", + "pd.DataFrame(post_tt_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"POST_ADM1_tt_ports.csv\"))\n", + "\n", + "post_zones_ports = ma.summarize_travel_time_populations(in_pop, post_friction, zones_df, mcp, tun_adm1)\n", + "pd.DataFrame(post_zones_ports.drop([\"geometry\"], axis=1)).to_csv(\n", + " os.path.join(results_folder, \"POST_ADM1_tt_zones.csv\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Zonal stats on nighttimelights" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "c:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\urllib3\\connectionpool.py:1099: InsecureRequestWarning: Unverified HTTPS request is being made to host 'globalnightlight.s3.amazonaws.com'. Adding certificate verification is strongly advised. See: https://urllib3.readthedocs.io/en/latest/advanced-usage.html#tls-warnings\n", + " warnings.warn(\n" + ] + } + ], + "source": [ + "ntl_files = ntlMisc.aws_search_ntl()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Clip out NTL raster files for pre and post\n", + "pre_file = [x for x in ntl_files if \"201501\" in x][0]\n", + "post_file = [x for x in ntl_files if \"202401\" in x][0]\n", + "\n", + "ntl_out_folder = os.path.join(base_folder, \"DATA\", \"NTL_Rasters\")\n", + "with rasterio.Env(GDAL_HTTP_UNSAFESSL = 'YES') as env:\n", + " pre_res = rMisc.clipRaster(rasterio.open(pre_file), tun_adm1, os.path.join(ntl_out_folder, \"VIIRS_201501.tif\"))\n", + " post_res = rMisc.clipRaster(rasterio.open(post_file), tun_adm1, os.path.join(ntl_out_folder, \"VIIRS_202401.tif\"))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "ntl_diff = post_res[0] - pre_res[0]\n", + "with rasterio.open(os.path.join(ntl_out_folder, \"VIIRS_201501_202401.tif\"), 'w', **pre_res[1]) as outR:\n", + " outR.write(ntl_diff)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Measure distance to treated roads" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "tPrint(\"Start\")\n", + "for road_file, out_file, col_lbl in [\n", + " [treated_segments_file, \"distance_to_treated_roads.csv\", \"RTE_NOM\"],\n", + " [control_segments_file, \"distance_to_control_roads.csv\", \"RTE_NOM\"],\n", + " [road_segments_file, \"distance_to_WB_digitized_roads.csv\", \"road_group\"],\n", + " ]:\n", + " roads = gpd.read_file(road_file)\n", + " roads = roads.to_crs(m_crs)\n", + " tun_adm2 = tun_adm2.to_crs(m_crs)\n", + " for lbl, df in roads.groupby(col_lbl):\n", + " tun_adm2[f'dist_{lbl}'] = tun_adm2.apply(lambda x: x[\"geometry\"].distance(df.union_all()), axis=1)\n", + " tun_adm2[f'dist_road'] = tun_adm2.apply(lambda x: x[\"geometry\"].distance(roads.union_all()), axis=1)\n", + " pd.DataFrame(tun_adm2.drop(['geometry'], axis=1)).to_csv(os.path.join(base_folder, \"RESULTS\", out_file))\n", + " tPrint(out_file)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gpd.read_file(control_segments_file)['RTE_NOM'].unique()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PRepare mapping data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "results_csv = os.path.join(base_folder, \"DATA\", \"MAPPING\", \"mapping_res.csv\")\n", + "in_res = pd.read_csv(results_csv)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "pd.merge(tun_adm1, in_res, left_on=\"shapeName\", right_on='adm1', how='outer').to_file( \n", + " os.path.join(base_folder, \"DATA\", \"MAPPING\", \"adm1_res.gpkg\"), driver=\"GPKG\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pd.merge(tun_adm1, in_res, left_on=\"shapeName\", right_on='adm1', how='outer').sort_values('time_hub')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gn", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..49f7c01 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,9 @@ +rasterio +geopandas +pandas +numpy +osmnx>1.3.0 +scikit-image +pyproj +scipy +pandana \ No newline at end of file