Create_accessibility_layers.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 creates accessibility surface and catchments for each point of interest in a PoI layer.\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 ...\n",
    "* **Population raster layer**: It should be in the same resolution with friction layers (100m). Building foorprint constrained population layers that WorldPop produces are suggested. They can be downloaded here:https://www.worldpop.org/project/categories?id=3\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",
    "* **Admin boundary**:  Geographic extent of area. It can be a country extent , a province or a district. It will be used to determine outher boundary of catchments. In other words, the catchments will be nested into specified  admin boundaries. \n",
    "\n",
    "## Output layers: \n",
    "\n",
    "Output layers will be saved into a directory that is created specified workspace directory.Two gdb will be created in the directory to save output layers\n",
    " <br>\n",
    "* **output.gbd** : The scripts run trought each specified admin units and creates output layers.It loops trought admin unit to nest output layers into admin units. Each admin units will have these layers ( \"...\" is place holder for name conventions that area specified in the script):<br>\n",
    "   ... _accessibility_... : Accessibility raster  layer in minute per a meter.<br>\n",
    "   ..._catchments_... :Catchments boundary<br>\n",
    "   ..._Sett_extent_cost_in_min_... : Settlement extent with distance to the closest PoI in minutes<br>\n",
    "   ..._sett_path_... : The shorthest path from a settlement to the closest PoI<br>\n",
    "   ..._sett_path_smoothlines_... : Smoothed version of the shorthest path<br>\n",
    " <br>\n",
    "* **final.gbd**: All indivual admin specific admin outputs will be merged and saved in this gdb.These layers will be created in the final.gdb:<br>\n",
    "   ..._based_accessibility_in_min: Accessibility raster  layer in minute per a meter.<br>\n",
    "   ..._based_catchments :Catchments boundary<br>\n",
    "   ...\"_based_settlements_in_min : Settlement extent with distance to the closest PoI in minutes<br>\n",
    "   ...\"_based_settlement_shortest_path : Smoothed version of the shorthest path<br>\n",
    "   <br>\n",
    "Note: Catchments can be created at a county extent ( without nesting into admin units). When this is the case, \"final.gdb\" will be empty. Output layers in the \"output.gdb\" should be used.    \n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "import numpy as np\n",
    "from arcgis.features import SpatialDataFrame\n",
    "arcpy.env.overwriteOutput = True\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": [
    "## User Inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "###=================== Initialize Variables ===================###\n",
    "\n",
    "country_iso=\"\"       # three letter country iso code e.g. zmb for Zombia\n",
    "poi_type=  \"\"           # type of point of interest such as \"health_facility\"      !! no space between names !!\n",
    "\n",
    "travel_type= \"\"       # options: walk or mix( walk+motorized)\n",
    "poi_subtype=  \"\"      # subtype of PoI  type. If there is no subtype , use \"all\"\n",
    "admin_col=\"\"         # field/column that hold admin names in admin boundary layer.\n",
    "                          # if you want to create output layers without nesting them into\n",
    "                            # a admin unit, create a dummy variable such as country iso code \n",
    "                            # then specify admin_col variale with the dummy variable\n",
    "\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 population raster\n",
    "path_to_pop_raster=r\"\" # WorldPop building footprint contraint population layer\n",
    "\n",
    "#path to admin  boundary\n",
    "path_to_pop_admin_boundary=r\"\" \n",
    "\n",
    "# path to settlement extent\n",
    "path_to_sett_extent=r\"\"\n",
    "### important: Settlement extent layers (bua, saa , hamlets) need to be merged. \n",
    "\n",
    "# path to the directory that you want to save outputs\n",
    "output_path=r\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "###==================================== functions =============================================###\n",
    "\n",
    "def project_raster(raster_to_be_projected, reference_raster):\n",
    "    '''project a raster layer'''\n",
    "    arcpy.env.snapRaster  =reference_raster\n",
    "    arcpy.env.extent      =reference_raster\n",
    "    arcpy.env.cellAlignment = \"ALIGN_WITH_PROCESSING_EXTENT\"\n",
    "    arcpy.env.cellSize = reference_raster\n",
    "    arcpy.env.overwriteOutput = True \n",
    "    arcpy.env.outputCoordinateSystem = arcpy.Describe(reference_raster).spatialReference\n",
    "    prj= arcpy.Describe(reference_raster).spatialReference\n",
    "    arcpy.ProjectRaster_management(raster_to_be_projected,\"pop_proj\",prj)\n",
    "    arcpy.env.workspace = output_gdb\n",
    "\n",
    "\n",
    "def spatial_join(target, join, fieldname, out_name):\n",
    "    '''spatial join two layers by taking mean of target variable'''\n",
    "    fieldmappings = arcpy.FieldMappings()\n",
    "    fieldmappings.addTable(target)\n",
    "    fieldmappings.addTable(join)\n",
    "    fieldIndex_ = fieldmappings.findFieldMapIndex(\"grid_code\")\n",
    "    fieldmap = fieldmappings.getFieldMap(fieldIndex_ )\n",
    "    field = fieldmap.outputField\n",
    "    field.name = fieldname\n",
    "    field.aliasName = fieldname\n",
    "    fieldmap.outputField = field\n",
    "    fieldmap.mergeRule = \"mean\"\n",
    "    fieldmappings.replaceFieldMap(fieldIndex_ , fieldmap)\n",
    "    arcpy.SpatialJoin_analysis(target, join,out_name, \"#\", \"#\", fieldmappings)\n",
    "\n",
    "\n",
    "def updated_boundary_to_admin(boundaryToBeAdjusted, admin_boundary, out_name):\n",
    "    '''nest catchments into admins boundary'''\n",
    "    smooth_boundaries=os.path.join(output_gdb,\"smooth_boundaries\")\n",
    "    clip_by_admin=os.path.join(output_gdb,\"clip_by_admin\")\n",
    "    union_with_admin=os.path.join(output_gdb,\"union_by_admin\")\n",
    "    multipart=os.path.join(output_gdb,\"multipart\")\n",
    "    neighbors_table=os.path.join(output_gdb,\"neighbors_table\")\n",
    "    neighbors_table_npArray=os.path.join(output_gdb,\"neighbors_table_npArray\")\n",
    "    neighbors_table_df=os.path.join(output_gdb,\"neighbors_table_df\")\n",
    "    multipart_npArray=os.path.join(output_gdb,\"multipart_npArray\")\n",
    "    multipart_df=os.path.join(output_gdb,\"multipart_df\")\n",
    "    multipart_temp=os.path.join(output_gdb,\"multipart_temp\")\n",
    "    multipart_final=os.path.join(output_gdb,out_name)\n",
    "    \n",
    "    ##======================================================================##\n",
    "    arcpy.cartography.SmoothPolygon(boundaryToBeAdjusted,smooth_boundaries,\n",
    "                                    \"PAEK\", \"200 Meters\", \"FIXED_ENDPOINT\", \"RESOLVE_ERRORS\", None)\n",
    "    \n",
    "    arcpy.Clip_analysis(smooth_boundaries,admin_boundary,clip_by_admin)\n",
    "    arcpy.Union_analysis([clip_by_admin,admin_boundary],union_with_admin)\n",
    "    arcpy.MultipartToSinglepart_management(union_with_admin,multipart)\n",
    "    \n",
    "    arcpy.AddField_management(multipart,\"unique_id\", \"TEXT\")\n",
    "    arcpy.CalculateField_management(multipart,\"uniqueid\",'str(!OBJECTID!)+\"_\"+str(!gridcode!)', \"PYTHON_9.3\",\"\")\n",
    "    \n",
    "    arcpy.PolygonNeighbors_analysis(multipart,neighbors_table, [\"uniqueid\"],\"\",\n",
    "                                       \"BOTH_SIDES\",\"\",\"METERS\",\"SQUARE_KILOMETERS\")\n",
    "    neighbors_table_npArray =arcpy.da.TableToNumPyArray(neighbors_table,\"*\")\n",
    "    neighbors_table_df=pd.DataFrame(neighbors_table_npArray)\n",
    "    \n",
    "    arcpy.AddField_management(multipart,\"area_km\", \"DOUBLE\")\n",
    "    arcpy.CalculateGeometryAttributes_management(multipart, [ [\"area_km\", \"AREA_GEODESIC\"]],\n",
    "                                                 \"\",\"SQUARE_KILOMETERS\")\n",
    "    \n",
    "    multipart_npArray=arcpy.da.FeatureClassToNumPyArray( multipart,[\"uniqueid\",\"area_km\"])\n",
    "    multipart_df=pd.DataFrame(multipart_npArray)\n",
    "    merge_tables=pd.merge(neighbors_table_df,multipart_df,left_on=\"src_uniqueid\", right_on=\"uniqueid\")\n",
    "    \n",
    "    merge_tables=merge_tables[merge_tables[\"area_km\"]<=0.1]\n",
    "   \n",
    "    merge_tables[\"gridcode2\"]=np.nan\n",
    "    get_unique_id=list(set(merge_tables.src_uniqueid.tolist()))\n",
    "\n",
    "    \n",
    "    for i in get_unique_id:\n",
    "        i_check=i.split(\"_\")[1:]\n",
    "       \n",
    "        rows=merge_tables.loc[merge_tables[\"src_uniqueid\"]==i]\n",
    "        rows_=rows[~rows[\"nbr_uniqueid\"].str.split(\"_\", expand=True)[1].isin(i_check)]\n",
    "        if rows_.shape[0]>=1:\n",
    "            get_max_lenght=rows_.LENGTH.idxmax()\n",
    "            get_max_lenght_class=rows.loc[get_max_lenght,\"nbr_uniqueid\"]\n",
    "            row_index=list(rows.index.values)\n",
    "            merge_tables.loc[row_index,'gridcode2_']=get_max_lenght_class\n",
    "         \n",
    "    merge_tables=merge_tables[~merge_tables[\"gridcode2_\"].isnull()] \n",
    "    merge_tables.to_csv(output_loc+\"\\\\merge_table.csv\") \n",
    "    merge_tables[[\"var1\",\"gridcode2\"]] = merge_tables['gridcode2_'].str.split(\"_\", expand=True)\n",
    "    merge_tables=merge_tables[[\"src_uniqueid\",\"nbr_uniqueid\",'gridcode2']]\n",
    "    sdf=pd.DataFrame.spatial.from_featureclass(multipart)[[\"gridcode\",\"uniqueid\",\"SHAPE\"]]\n",
    "    sdf_merge=pd.merge(sdf,merge_tables,left_on=\"uniqueid\",right_on=\"src_uniqueid\",how=\"left\")\n",
    "    sdf_merge.loc[sdf_merge['gridcode2'].isnull(),'gridcode2'] = sdf_merge['gridcode']   \n",
    "    sdf_merge.spatial.to_featureclass(multipart_temp)\n",
    "    arcpy.Dissolve_management(multipart_temp,multipart_final,[[\"gridcode2\",\"uniqueid\"]],\n",
    "                              None, \"SINGLE_PART\", \"DISSOLVE_LINES\") \n",
    "    arcpy.AlterField_management(multipart_final, \"gridcode2\", 'gridcode',\"gridcode\")\n",
    "    ##===========================================================================###\n",
    "    \n",
    "    delete_fc=[smooth_boundaries,clip_by_admin,union_with_admin,multipart,multipart_temp]\n",
    "  \n",
    "    for fc in delete_fc:\n",
    "        arcpy.Delete_management(fc) \n",
    "    \n",
    "\n",
    "\n",
    "def create_accessibility_layers2(admin_name):\n",
    "    ''' creates catchments for a PoI layer''' \n",
    "    \n",
    "    if not arcpy.Exists(admin_name+\"_accessibility_\"+poi_subtype):\n",
    "        \n",
    "        ## initiliaze admin specific temporal layers\n",
    "        admin_boundary=admin_name+\"_boundary\"\n",
    "        admin_poi=admin_name+\"_poi\"\n",
    "        admin_sett=admin_name+\"_sett\"\n",
    "        admin_friction=admin_name+\"_friction\"\n",
    "        \n",
    "        ##==========================================================================================###\n",
    "        ## initialize output layer names\n",
    "        cost_distance=admin_name+\"_accessibility_\"+  poi_subtype\n",
    "        catchment_areas_update= admin_name+\"_catchments_\"+  poi_subtype\n",
    "        sett_path=admin_name+\"_sett_path_\"+  poi_subtype\n",
    "        sett_path_smoothline=admin_name+\"_sett_path_smoothlines_\"+  poi_subtype\n",
    "        sett_with_cost=admin_name+\"_Sett_extent_cost_in_min_\"+  poi_subtype\n",
    "\n",
    "        ##==========================================================================================### \n",
    "\n",
    "        ## extract layers by the admin unit\n",
    "        arcpy.Select_analysis(admins_bndry,admin_boundary,\"admin_cln = '{}'\".format(admin_name))\n",
    "        arcpy.Select_analysis(poi,admin_poi,\"admin_cln = '{}'\".format(admin_name))\n",
    "        arcpy.gp.ExtractByMask_sa(friction_layer,admin_boundary,admin_friction)\n",
    "        arcpy.Clip_analysis(path_to_sett_extent, admin_boundary,admin_sett)\n",
    "        \n",
    "        if int(arcpy.GetCount_management(admin_poi).getOutput(0))>=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",
    "            print ( \"   >>> Creating catcment areas of PoI  based on  minimum cost\")\n",
    "            arcpy.sa.CostAllocation(admin_poi,admin_friction,\"\",\"\",\"OBJECTID\",\"cost_allocation\").save(\"cost_allocation_temp\")\n",
    "            arcpy.RasterToPolygon_conversion(\"cost_allocation_temp\", \"catchment_areas_temp\", \"NO_SIMPLIFY\")\n",
    "\n",
    "            #------------------------------------------------------------------------------------------#    \n",
    "            print (\"   >>> Calculating population of each catchmet area by 30 min dist. interval\")     \n",
    "            arcpy.sa.Con(cost_distance,130, cost_distance, \"VALUE>120\").save(cost_distance+\"_con\")\n",
    "            arcpy.sa.Contour(cost_distance+\"_con\", cost_distance+\"_contour\", 30, 0, 1, \"CONTOUR_POLYGON\")\n",
    "            arcpy.AddField_management(cost_distance+\"_contour\", \"time_interval\", \"TEXT\")\n",
    "            arcpy.management.CalculateField(cost_distance+\"_contour\", \"time_interval\", 'str(!ContourMin!)+\"-\"+str(!ContourMax!)', \"PYTHON3\", '', \"TEXT\")\n",
    "            arcpy.Intersect_analysis([cost_distance+\"_contour\",\"catchment_areas_temp\"], \"catchment_areas_int\", \"ALL\")\n",
    "            arcpy.AddField_management( \"catchment_areas_int\", \"zone_unique\", \"TEXT\")\n",
    "            arcpy.management.CalculateField( \"catchment_areas_int\", \"zone_unique\", 'str(!gridcode!)+\"_\"+str(!time_interval!)', \"PYTHON3\", '', \"TEXT\")\n",
    "            arcpy.sa.ZonalStatisticsAsTable(\"catchment_areas_int\", \"zone_unique\", path_to_pop_raster, \"zonal_pop_min\",\"DATA\", \"SUM\")\n",
    "            totpop_min_array = arcpy.da.TableToNumPyArray(\"zonal_pop_min\",[\"zone_unique\",\"SUM\"])\n",
    "            df=pd.DataFrame(totpop_min_array)\n",
    "            df[['gridcode','time_interval']] = df[ \"zone_unique\"].str.split(\"_\",expand=True)\n",
    "            df= df.pivot(index='gridcode', columns='time_interval', values=\"SUM\")\n",
    "            df.rename({\"0.0-30.0\":\"under_30min\",\"30.0-60.0\":\"30min_60min\",\"60.0-90.0\":\"60min_90min\",\"90.0-120.0\":\"90min_120min\",\n",
    "                         \"120.0-130.0\":\"above_120min\"}, axis=1, inplace=True)\n",
    "            df=   df.fillna(0)\n",
    "            col_to_sum=[c for c in df.columns if  c.endswith(\"min\")]\n",
    "            df[\"totpop_20\"]=df[col_to_sum].sum(axis=1)\n",
    "            df_=pd.DataFrame(df.to_records())\n",
    "            df_[\"gridcode\"]=df_[\"gridcode\"].astype(str).astype(\"int64\")\n",
    "           \n",
    "\n",
    "            print (\"   >>> Updating catchment boundary based on admin boundary\")\n",
    "            updated_boundary_to_admin(\"catchment_areas_temp\", admin_boundary,catchment_areas_update)\n",
    "            \n",
    "            sdf=pd.DataFrame.spatial.from_featureclass(catchment_areas_update)\n",
    "            sdf[\"gridcode\"]=sdf[\"gridcode\"].astype(str).astype(\"int64\")\n",
    "            sdf_poi=pd.DataFrame.spatial.from_featureclass(admin_poi)\n",
    "            sdf_poi.drop([\"SHAPE\", 'NEAR_DIST',\"NEAR_FID\",\"admin_cln\"], axis=1, inplace=True)\n",
    "            sdf_poi.rename({\"OBJECTID\":\"gridcode\"},axis=1, inplace=True)\n",
    "       \n",
    "            \n",
    "            sdf_merge1=sdf.merge(df_, how=\"left\",on=\"gridcode\")\n",
    "            sdf_merge2=sdf_merge1.merge(sdf_poi, how=\"left\",on=\"gridcode\")\n",
    "            numeric_columns = sdf_merge2.select_dtypes(include=['int64',\"float64\"]).columns\n",
    "            sdf_merge2[numeric_columns ]=sdf_merge2[numeric_columns].fillna(0)\n",
    "            numeric_columns = sdf_merge2.select_dtypes(include=['object']).columns\n",
    "            sdf_merge2[numeric_columns ]=sdf_merge2[numeric_columns].fillna('')\n",
    "            arcpy.Delete_management(catchment_areas_update)\n",
    "            sdf_merge2.spatial.to_featureclass(output_gdb+\"\\\\\"+catchment_areas_update)\n",
    "\n",
    "            #---------------------------------------------------------------------------------#\n",
    "\n",
    "            print (\"   >>> Creating cost path from settlements to PoI\")\n",
    "            arcpy.sa.CostPathAsPolyline(path_to_sett_extent, cost_distance, \"outBkLinkRaster\", sett_path, \"EACH_ZONE\", \"OBJECTID\", \"INPUT_RANGE\")\n",
    "\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",
    "            arcpy.DeleteField_management( sett_with_cost,  [\"Join_Count\", \"TARGET_FID\", \"pointid\"])\n",
    "\n",
    "            print (\"   >>> Deleting temporary layers...\")\n",
    "            delete_fc=[\"outBkLinkRaster\",\"cost_allocation_temp\", \"cost_distance_point\",\"catchment_areas_int\", \"catchment_areas_temp\",\"zonal_pop_min\",\n",
    "                       \"settlement_path_temp\",\"cost_allocation\",\"cost_distance_all_con\",\"cost_distance_all_contour\",admin_name+\"_boundary\",\n",
    "                       admin_name+\"_poi\",admin_name+\"_poi\",admin_name+\"_friction\",admin_name+\"_sett\",admin_name+\"_accessibility_all_con\",\n",
    "                      admin_name+\"_accessibility_all_contour\" ]\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 admin specif output vector layers'''\n",
    "    merge_fc=[fc for fc in arcpy.ListFeatureClasses() if fc.endswith(select_layer_to_be_merged)]\n",
    "    if len(merge_fc):\n",
    "        arcpy.Merge_management(merge_fc, final_gdb+\"\\\\\"+output_name)\n",
    "\n",
    "def merge_raster(select_layer_to_be_merged, output_name):\n",
    "    '''merge admin specif output raster layers'''\n",
    "    merge_rasters=[fc for fc in arcpy.ListRasters() if fc.endswith(\"cost_distance_all\")]\n",
    "    if len(merge_rasters):\n",
    "        arcpy.MosaicToNewRaster_management(merge_rasters, final_gdb, output_name,\"\",\"\",\"\",1)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "###================================ create output workspace ==========================================###\n",
    "\n",
    "# create output directory\n",
    "if not  os.path.exists (os.path.join(output_path,\"output_\"+poi_type+\"_\"+poi_subtype+\"_\"+travel_type)):\n",
    "    os.mkdir(os.path.join(output_path,\"output_\"+poi_type+\"_\"+poi_subtype+\"_\"+travel_type))\n",
    "output_loc=os.path.join(output_path,\"output_\"+poi_type+\"_\"+poi_subtype+\"_\"+travel_type)\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",
    "\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",
    "arcpy.env.workspace = output_gdb\n",
    "\n",
    "###======================= copy admin units and PoI into workspace ==================================###\n",
    "\n",
    "arcpy.CopyFeatures_management(path_to_poi, poi_type+\"_\"+poi_subtype)\n",
    "poi=poi_type+\"_\"+poi_subtype\n",
    "arcpy.CopyFeatures_management(path_to_pop_admin_boundary, \"admin_boundary\")\n",
    "admins_bndry=\"admin_boundary\"\n",
    "poi=\"hf_all\"\n",
    "\n",
    "\n",
    "\n",
    "#==================================== format admin names =============================================###\n",
    "objectid_list=[]\n",
    "admin_list_=[]\n",
    "arcpy.AddField_management(admins_bndry, \"admin_cln\", \"TEXT\")\n",
    "with arcpy.da.UpdateCursor(admins_bndry,[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",
    "\n",
    "admin_list2_=[]\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(poi,admins_bndry)\n",
    "arcpy.AddField_management(poi, \"admin_cln\", \"TEXT\")\n",
    "for i, name_ in zip(objectid_list,admin_list_): \n",
    "    with arcpy.da.UpdateCursor(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(name_)\n",
    "            cursor.updateRow(row)\n",
    "\n",
    "admin_list=list(set(admin_list2_))\n",
    "\n",
    "#==================================== project population raster ==================================###\n",
    "\n",
    "project_raster(path_to_pop_raster, friction_layer)\n",
    "path_to_pop_raster=\"pop_proj\"\n",
    "\n",
    "\n",
    "#==================================== main process =============================================###\n",
    "\n",
    "\n",
    "\n",
    "print(f\"There are {len(admin_list)} admin units to be processed...\" )\n",
    "\n",
    "## Create output by each admin units\n",
    "for admin in admin_list:\n",
    "\n",
    "    try:\n",
    "        print (\"##########  \"+admin+\"  ###########\")\n",
    "        create_accessibility_layers2(admin)\n",
    "    except:\n",
    "        print (admin+\"  did NOT RUN\")\n",
    "\n",
    "#=========================== merge admin specific admin outputs ==========================###\n",
    "(\"   >>> Merging district based output layers...\")\n",
    "merge_outputs(\"catchment_update_all\",country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_based_catchments\" )\n",
    "merge_outputs(\"sett_path_smoothlines_all\",country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_based_settlement_shortest_path\" )\n",
    "merge_outputs(\"Sett_extent_cost_in_min_all\",country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_based_settlements_with_time\" )\n",
    "merge_raster(\"cost_distance_all\", country_iso+\"_\"+poi_type+\"_\"+travel_type+\"_based_access_surface_in_min\")\n",
    "\n",
    "print (\"SCRIPT IS DONE\")\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "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
}