Calculate_accessibility_by_distace.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Hasim Engin <br>\n",
    "email:hengin@ciesin.columbia.edu<br>\n",
    "Project: GRID3<br>\n",
    "Organization:CIESIN, Columbia University<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create Catchments for a PoI "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook calculates dictance from s settlement to the cloesest PoI by linear distance based on a friction surface.\n",
    "## Input layers:\n",
    "\n",
    "* **Friction layers**: Friction layers (walk and walk+motorized) in minutes.They can be created by running \"create_friction_surface\" script.\n",
    "* **PoI layer**: A point of interest layer. Such as health facilities, schools, water sources .... POI layer must have unique id and capacity for each record\n",
    "* **Settlement extent layer**: This is a vector layer produced by GRID3 for sub-saharas countries. They can be downloaded here:https://grid3.org/resources/data\n",
    "\n",
    "\n",
    "\n",
    "## Output layers: \n",
    "\n",
    "Output layers will be saved into a directory that is created specified workspace directory.The scripts goes trought individual admin units that is specified and calculates distance between each individual settlement to the closest PoI.\n",
    "Individual admin units result will be saved in the output.gdb. Then all the admin level settlement extent will merged and saved in the final.gdb. Output.gdb can be deleted if output layer in the final.gdb is good.\n",
    " <br>\n",
    "* **output.gbd** : admin level settlement extent will be saved here. \n",
    " <br>\n",
    "* **Final.gdb** : <br>\n",
    "   ..._settlement_extent: Settlement extent with distance to closest PoI in meter.<br>\n",
    "   ..._settlement_part: The shortest path that connects a settlement to a health facility. Distance is calculated based on these paths. \n",
    "   \n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# import libraries\n",
    "import os\n",
    "import arcpy\n",
    "import re\n",
    "import unidecode\n",
    "import pandas as pd\n",
    "import arcpy.cartography as CA\n",
    "from arcgis.features import SpatialDataFrame\n",
    "arcpy.env.overwriteOutput = True\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# required libraries can be installed with following code\n",
    "#uncommaent the line below and write a -library name- after !):\n",
    "#! -library name- install"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "###==================================== functions =============================================###\n",
    "\n",
    "def calculate_shorthest_distrance(admin_name):\n",
    "    admin_name=str(admin_name)\n",
    "        ##Creates cathment areas along with some popualtion summary\n",
    "        ##Cathments areas are nested in the admin units\n",
    "       \n",
    "    if not arcpy.Exists(\"settlement_extent_\"+admin_name):\n",
    "       \n",
    "        admin_boundary=\"boundary_\"+admin_name\n",
    "        admin_poi=\"poi_\"+admin_name\n",
    "        admin_sett=\"sett_\"+admin_name\n",
    "        admin_sett_point=\"sett_point_\"+admin_name\n",
    "        admin_friction=\"friction_\"+admin_name\n",
    "        \n",
    "        \n",
    "        ##==========================================================================================###\n",
    "        cost_distance=\"accessibility_\"+admin_name\n",
    "        catchment_areas_update=\"catchments_\"+  admin_name+\"_\"+poi_subtype\n",
    "        sett_path=\"sett_path_\"+ admin_name+\"_\"+poi_subtype\n",
    "        sett_path_smoothline=\"sett_path_smoothlines_\"+admin_name+\"_\"+poi_subtype\n",
    "        sett_with_cost=\"Sett_extent_cost_in_min_\"+ admin_name +\"_\"+poi_subtype\n",
    "\n",
    "        ##==========================================================================================### \n",
    "\n",
    "        ## extract layers by the admin unit\n",
    "        arcpy.Select_analysis(path_to_admin_boundary,admin_boundary,\"OBJECTID = {}\".format(admin_name))\n",
    "        arcpy.analysis.Buffer(admin_boundary,\"buffer\", \"100 Meters\", \"FULL\", \"ROUND\", \"ALL\", None, \"PLANAR\")\n",
    "        arcpy.Select_analysis(path_to_poi,admin_poi,\"NEAR_FID = {}\".format(admin_name))\n",
    "\n",
    "        arcpy.gp.ExtractByMask_sa(friction_layer,\"buffer\",admin_friction)\n",
    " \n",
    "        \n",
    "        arcpy.Clip_analysis(path_to_sett_extent, admin_boundary,admin_sett)\n",
    "        arcpy.FeatureToPoint_management(admin_sett,admin_sett_point, \"INSIDE\")\n",
    "\n",
    "        check_bdry=int(arcpy.GetCount_management(admin_boundary).getOutput(0))\n",
    "        point_check=int(arcpy.GetCount_management(admin_poi).getOutput(0))\n",
    "        print (check_bdry,point_check)\n",
    "        if  check_bdry>=1 and point_check>=1 :\n",
    "\n",
    "\n",
    "            print (f\"Creating accessiblility layers for =={country_iso}== for =={poi_type}== with =={travel_type}== travel mode\" )\n",
    "            print(\" \")\n",
    "            print ( \"   >>> Creating cost distance in minutes to PoI\")\n",
    "            arcpy.sa.CostDistance(admin_poi,admin_friction,\"\" , \"outBkLinkRaster\").save(  cost_distance)\n",
    "\n",
    "\n",
    "            #---------------------------------------------------------------------------------#\n",
    "\n",
    "            print (\"   >>> Creating cost path from settlements to PoI\")\n",
    "            arcpy.sa.CostPathAsPolyline(admin_sett_point, cost_distance, \"outBkLinkRaster\", sett_path, \"EACH_ZONE\", \"OBJECTID\", \"INPUT_RANGE\")\n",
    "            arcpy.AddField_management(sett_path, \"dist_in_meter\", \"DOUBLE\")\n",
    "            arcpy.CalculateGeometryAttributes_management(sett_path, \"dist_in_meter LENGTH_GEODESIC\", \"METERS\")\n",
    "            sdf_path=pd.DataFrame.spatial.from_featureclass(sett_path)[[\"PathCost\",\"DestID\",\"dist_in_meter\"]]\n",
    "            sdf_path.rename(columns={\"PothCost\":\"dist_in_min\",\"DestID\":\"matchid\"}, inplace=True)\n",
    "            print (\"   >>> Smoothing paths\")\n",
    "            CA.SimplifyLine(sett_path,sett_path_smoothline, \"POINT_REMOVE\", 100,\"\",\"NO_KEEP\")\n",
    "\n",
    "            print ( \"   >>> Calculating everage travel time in minutes to closest PoI...\")\n",
    "            arcpy.RasterToPoint_conversion(cost_distance, \"cost_distance_point\", \"VALUE\")\n",
    "            spatial_join(admin_sett, \"cost_distance_point\", \"average_cost_min\", sett_with_cost)\n",
    "            sdf_sett=pd.DataFrame.spatial.from_featureclass(sett_with_cost)[[\"OBJECTID_1\",\"OBJECTID\",\"MGRS_Code\",\"average_cost_min\", \"SHAPE\"]]\n",
    "            sdf_sett.rename(columns={\"OBJECTID\":\"matchid\",\"average_cost_min\":\"avrg_dist_min\"}, inplace=True)\n",
    "            sdf_merge=sdf_sett.merge(sdf_path, how=\"left\",on=\"matchid\")\n",
    "            sdf_merge.spatial.to_featureclass(output_gdb+\"\\\\\"+\"settlement_extent_\"+admin_name)\n",
    "\n",
    "            print (\"   >>> Deleting temporary layers...\")\n",
    "            delete_fc=[\"sett_point_\"+admin_name,\"boundary_\"+admin_name,\"buffer\",\"friction_\"+admin_name,\"outBkLinkRaster\",\"poi_\"+admin_name,\"sett_\"+admin_name,\n",
    "                       \"accessibility_\"+admin_name,\"Sett_extent_cost_in_min_\"+admin_name+\"_all\"]\n",
    "            for fc in delete_fc:\n",
    "                 arcpy.Delete_management(fc)      \n",
    "            print (\"##====================================================================================================##\")    \n",
    "\n",
    "\n",
    "def merge_outputs(select_layer_to_be_merged, output_name):\n",
    "    merge_fc=[fc for fc in arcpy.ListFeatureClasses() if fc.startswith(select_layer_to_be_merged)]\n",
    "    if len(merge_fc)>1:\n",
    "        arcpy.Merge_management(merge_fc, final_gdb+\"\\\\\"+output_name)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "###=================== Initialize Variables ===================###\n",
    "\n",
    "country_iso=\"\"          # three latters iso code\n",
    "poi_type=  \"\"           # type of point of interest such as \"health_facility\"      !! no space between names !!\n",
    "travel_type= \"\"         # options: walk or mix( walk+motorized)                      \n",
    "admin_col=\"\"            # admin column in admin boundary layer.\n",
    "###================================ input layers ==============================================###\n",
    "\n",
    "# path to PoI layer such as health facilities or schools\n",
    "path_to_poi=r\"\"\n",
    "\n",
    "# friction surface for walking. It is created by \"create_friction_surface\" script\n",
    "path_to_friction_walk=r\"\"\n",
    "\n",
    "# friction surface for motorized.It is created by \"create_friction_surface\" script\n",
    "path_to_friction_mix=r\"\"\n",
    "\n",
    "#path to settlement extent layer\n",
    "path_to_sett_extent==r\"\" \n",
    "\n",
    "#path to admin  boundary\n",
    "path_to_admin_boundary=r\"\"\n",
    "\n",
    "# path to the directory that you want to save outputs\n",
    "output_path=r\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Output Workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    " #create output directory\n",
    "if not  os.path.exists (os.path.join(output_path,\"output\")):\n",
    "    os.mkdir(os.path.join(output_path,\"output\"))\n",
    "output_loc=os.path.join(output_path,\"output\")\n",
    "\n",
    " #create output gdb\n",
    "if  not arcpy.Exists(os.path.join(output_loc,\"output.gdb\")):\n",
    "    arcpy.CreateFileGDB_management(output_loc,\"output.gdb\")\n",
    "output_gdb=os.path.join(output_loc,\"output.gdb\")\n",
    " #create final gdb\n",
    "if  not arcpy.Exists(os.path.join(output_loc,\"final.gdb\")):\n",
    "    arcpy.CreateFileGDB_management(output_loc,\"final.gdb\")\n",
    "final_gdb=os.path.join(output_loc,\"final.gdb\")\n",
    "\n",
    "# workspace\n",
    "arcpy.env.workspace=output_gdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "###==================== select friction layer based on travel type============================###\n",
    "\n",
    "if travel_type==\"walk\":\n",
    "    friction_layer= path_to_friction_walk\n",
    "if travel_type==\"mix\":\n",
    "    friction_layer= path_to_friction_mix\n",
    "    \n",
    "###=================== transfer admin name from admin boundary to PoI =======================###   \n",
    "# get admin names ( no special characters)   and respective abjectids\n",
    "objectid_list=[]\n",
    "admin_list=[]\n",
    "arcpy.AddField_management(path_to_admin_boundary, \"admin_cln\", \"TEXT\")\n",
    "with arcpy.da.UpdateCursor(path_to_admin_boundary,[admin_col,\"admin_cln\", 'OBJECTID']) as cursor:\n",
    "    for row in cursor:\n",
    "        if row[0] is not None:\n",
    "            fix_admin= unidecode.unidecode(row[0])\n",
    "            fix_admin=re.sub(r'[-|$|.|!|)|(|!|-|_|]',r'',fix_admin) # spatial characters need to be removed in admin names\n",
    "            fix_admin=fix_admin.strip().replace(\"  \",\"_\").replace(\" \",\"_\")\n",
    "            row[1]=fix_admin\n",
    "            admin_list_.append(row[1])\n",
    "            objectid_list.append(row[2])\n",
    "            cursor.updateRow(row)\n",
    "\n",
    "            \n",
    " ##Transfer admin names from admin boundary to PoI layer.\n",
    " ##Admin names in boudary laeyr and PoI layer has to be identical\n",
    "arcpy.Near_analysis(path_to_poi,path_to_admin_boundary,\"10 Centimeters\")\n",
    "arcpy.AddField_management(path_to_poi, \"admin_cln\", \"TEXT\")\n",
    "for i, name_ in zip(objectid_list,admin_list_): \n",
    "    with arcpy.da.UpdateCursor(path_to_poi,[\"NEAR_FID\",\"admin_cln\"]) as cursor:\n",
    "        for row in cursor:\n",
    "            if row[0]==i:\n",
    "                row[1]=name_\n",
    "                admin_list2_.append(i)\n",
    "            cursor.updateRow(row)\n",
    "\n",
    "admin_list=list(set(admin_list)) \n",
    "###======================== prosessing ==============================###\n",
    "\n",
    "not_run=[]\n",
    "\n",
    "count=len(admin_list)\n",
    "## Create output by each admin units\n",
    "for admin in admin_list:\n",
    "    print(count)\n",
    "    try:\n",
    "        print (\"##########  \"+str(admin)+\"  ###########\")\n",
    "        create_accessibility_layers2(admin)\n",
    "    except:\n",
    "        print (str(admin)+\"  did NOT RUN\")\n",
    "        not_run.append(admin)\n",
    "    count=count-1\n",
    "\n",
    "(\"   >>> Merging district based output layers...\")\n",
    "merge_outputs(\"sett_path_smoothlines_\",country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_settlement_path\" )\n",
    "merge_outputs(\"settlement_extent_\",country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_settlement_extent\" )\n",
    "\n",
    "\n",
    "# get if any admin unit did not run\n",
    "if len (not_run) >=1:\n",
    "    print (\"These admin units did not run properly\")\n",
    "    for i in not_run:\n",
    "        print (i)\n",
    "print (\"SCRIPT IS DONE\")\n",
    "    \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ArcGISPro",
   "language": "Python",
   "name": "python3"
  },
  "language_info": {
   "file_extension": ".py",
   "name": "python",
   "version": "3"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}