Copy of storeygard_replication_assignment.ipynb

{"cells":[{"cell_type":"markdown","id":"2d8ba2db-a76c-427d-9dbb-71afc152cbb6","metadata":{"id":"2d8ba2db-a76c-427d-9dbb-71afc152cbb6"},"source":["## -------------------------------------------------------------------------\n","## YOUR SUBMISSION WILL NOT BE RECORDED UNLESS YOU COMPLETE YOUR DETAILS HERE\n","### *replace the text in the brackets below with your own details*\n","\n","### NAME: [First name, Last name]\n","### EMAIL: [email address that is registered with Wheeler]\n","\n","## -------------------------------------------------------------------------"]},{"cell_type":"markdown","id":"17baa446-b106-4e3d-b987-7e0eedfca88e","metadata":{"id":"17baa446-b106-4e3d-b987-7e0eedfca88e"},"source":["## Simplyfied Replication of Storeygard (2016) \n","\n","In this assignment you have to generate the (simplyfied) final dataset used in **\"Farther on down the road: transport costs, trade and urban growth\"** (2016) Review of Economic Studies 83(3): 1263-1295. For this assignment, we are going to focus on **Tanzania** for the year **1992** and **2013**.\n","\n","The final dataset (6 columns, 148 raws) should contained:\n","   - 74 Tanzanian city names (repeated twice for 1992 and 2013)\n","   - city populations (same value for 1992 and 2013)\n","   - a column for the year (1992 and 2013)\n","   - light values in both years (1992 and 2013)\n","   - oil prices in both years (1992 and 2013) \n","   - distance to the primary city (same value for 1992 and 2013)\n","\n","The objective is to use the geoprocessing tools and routines seen in class to solve the intermediate steps needed to assemble the final dataset. \n","\n","We have organized the assigment into 6 blocks:\n","\n","1. Identify Tanzanian Cities\n","2. Compute total light emitted by each city in 1992 and 2013\n","3. Identify Tanzania's primary city\n","4. Create a routable road network\n","5. Find optimal paths between cities\n","6. Build the final dataset\n","\n","Each section and subsection provide a detailed explanation of\n","\n","- Inputs: shapefiles, raster\n","- Task objective and description\n","- Output\n","\n","We are providing the packages we exepect you to use as well as the Coordinate Reference Systems.\n","\n","Good luck!\n"]},{"cell_type":"code","source":["%pip install geopandas"],"metadata":{"id":"wtdPVnYQx_g3","executionInfo":{"status":"ok","timestamp":1685888494664,"user_tz":-180,"elapsed":7631,"user":{"displayName":"Martin Sure Kondiwa","userId":"17309152213992783019"}},"outputId":"445b77b3-62a0-4411-ae56-1c83b531b772","colab":{"base_uri":"https://localhost:8080/"}},"id":"wtdPVnYQx_g3","execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n","Requirement already satisfied: geopandas in /usr/local/lib/python3.10/dist-packages (0.13.0)\n","Requirement already satisfied: fiona>=1.8.19 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.9.4.post1)\n","Requirement already satisfied: packaging in /usr/local/lib/python3.10/dist-packages (from geopandas) (23.1)\n","Requirement already satisfied: pandas>=1.1.0 in /usr/local/lib/python3.10/dist-packages (from geopandas) (1.5.3)\n","Requirement already satisfied: pyproj>=3.0.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (3.5.0)\n","Requirement already satisfied: shapely>=1.7.1 in /usr/local/lib/python3.10/dist-packages (from geopandas) (2.0.1)\n","Requirement already satisfied: attrs>=19.2.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (23.1.0)\n","Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (2022.12.7)\n","Requirement already satisfied: click~=8.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (8.1.3)\n","Requirement already satisfied: click-plugins>=1.0 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.1.1)\n","Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (0.7.2)\n","Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from fiona>=1.8.19->geopandas) (1.16.0)\n","Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2.8.2)\n","Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (2022.7.1)\n","Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.1.0->geopandas) (1.22.4)\n"]}]},{"cell_type":"code","execution_count":null,"id":"fa25d22b-eac0-4158-b052-a7e72bc9ae70","metadata":{"tags":[],"id":"fa25d22b-eac0-4158-b052-a7e72bc9ae70"},"outputs":[],"source":["import numpy as np\n","import pandas as pd\n","import geopandas as gpd\n","from pyproj import CRS\n","import rasterio\n","from rasterstats import zonal_stats\n","import matplotlib.pyplot as plt\n","from shapely.ops import nearest_points\n","from shapely.geometry import LineString, Point\n","import networkx as nx"]},{"cell_type":"code","source":["from google.colab import drive\n","drive.mount(\"/content/drive\")"],"metadata":{"id":"ZSmkzm0OydBY","executionInfo":{"status":"ok","timestamp":1685888626379,"user_tz":-180,"elapsed":3752,"user":{"displayName":"Martin Sure Kondiwa","userId":"17309152213992783019"}},"outputId":"ada7db7a-fec9-487a-d184-aac7aa2b5377","colab":{"base_uri":"https://localhost:8080/"}},"id":"ZSmkzm0OydBY","execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"]}]},{"cell_type":"code","source":["World=gpd.read_file('/content/drive/MyDrive/GIS_Research_2023/2/data/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp')"],"metadata":{"id":"ONq2XOS4zNoM"},"id":"ONq2XOS4zNoM","execution_count":null,"outputs":[]},{"cell_type":"code","source":["dir_path='/content/drive/MyDrive/GIS_Research_2023/Assignment'"],"metadata":{"id":"-g7vnW0q2nkX"},"id":"-g7vnW0q2nkX","execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"id":"2bf0910b-1ffd-4c51-a650-77c73557beac","metadata":{"id":"2bf0910b-1ffd-4c51-a650-77c73557beac"},"outputs":[],"source":["# Defining Coordinate Reference System to use throughout the whole notebook\n","wgs84_crs = CRS.from_string('EPSG:4326') # WGS 1984\n","aeaa_crs = CRS.from_string(\"esri:102022\") # albers equal area africa crs\n","aedc_crs = CRS.from_string(\"esri:102023\") # Africa equidistant conic"]},{"cell_type":"markdown","id":"4d6313bd-182f-4d05-9548-8fa7a6e49e9e","metadata":{"id":"4d6313bd-182f-4d05-9548-8fa7a6e49e9e"},"source":["# 1. Identify Tanzanian cities\n","\n","The unit of observations is the city. Later on, we would to like to calculate average luminosity in both 1992 and 2013 for each city. To do that we need to have a polygon associated with each city. We need to find a way to assign a city name to its corresponding polygon.\n","\n","- We have provided point locations of cities around the world (`citypop_v4_latlons.csv`) and their maximal night lights extents (`city_extents.shp`)\n","\n","- Write code to identify, in Tanzania, light extents with known cities and discard light extends that cannot be identified with a city\n","\n","\n"]},{"cell_type":"code","execution_count":null,"id":"da4fb13f-4835-4711-88b1-8e7d39933446","metadata":{"id":"da4fb13f-4835-4711-88b1-8e7d39933446"},"outputs":[],"source":["#import necessary libraries \n","import geopandas as gpd\n","import pandas as pd\n"]},{"cell_type":"code","source":["#Load city locations data\n","df=pd.read_csv(dir_path+'/data/citypop_v4_latlons.csv')"],"metadata":{"id":"Xha_peDc3S2M"},"id":"Xha_peDc3S2M","execution_count":null,"outputs":[]},{"cell_type":"code","source":["#Load city extents shapefile\n","city_extents_gdf=gpd.read_file(city_extents.shp)"],"metadata":{"id":"xWwqeMMQ5Ydg"},"id":"xWwqeMMQ5Ydg","execution_count":null,"outputs":[]},{"cell_type":"code","source":["#Filter extents to Tanzania only\n","tanzania_extents_gdf=city_extents_gdf[city_extents_gdf['COUNTRY']=='Tanzania']"],"metadata":{"id":"bOp6qMKa7qzp"},"id":"bOp6qMKa7qzp","execution_count":null,"outputs":[]},{"cell_type":"code","source":["#Join city locations data with extents on city name\n","tanzania_extents_gdf=tanzania_extents_gdf.join(cities_df.set_index('City'),on='city_name')"],"metadata":{"id":"GuX9PvMn8mqy"},"id":"GuX9PvMn8mqy","execution_count":null,"outputs":[]},{"cell_type":"markdown","source":[],"metadata":{"id":"jsNFVhJ4-OmO"},"id":"jsNFVhJ4-OmO"},{"cell_type":"code","source":["#Drop extents with no matching city\n","tanzania_cities_extents_gdf.dropna(subset=['City'],inplace=True)"],"metadata":{"id":"1Baay-1O9ST9"},"id":"1Baay-1O9ST9","execution_count":null,"outputs":[]},{"cell_type":"code","source":["#Export filtered extents to new shapefile"],"metadata":{"id":"5DgPht_495As"},"id":"5DgPht_495As","execution_count":null,"outputs":[]},{"cell_type":"code","source":["tanzania_extents_gdf.to_file('tanzania_city_extents.shp')"],"metadata":{"id":"YwByJ1YV-CO1"},"id":"YwByJ1YV-CO1","execution_count":null,"outputs":[]},{"cell_type":"markdown","id":"3028979b-c802-4a04-835d-e6d2f66a207d","metadata":{"id":"3028979b-c802-4a04-835d-e6d2f66a207d"},"source":["# 2. Compute total light emitted by each city in 1992 and 2013\n","\n","Now that you manage to assign each city point to a polygon, your task is to calculate the average luminosity level for each polygon for the years 1992 and 2013.\n","\n","- We have provided raster data of nighttime lights in Tanzania for 1992 and 2013 (`tza_lights_1992.tif` and `tza_lights_2013.tif`). \n","\n","- Combine them with the identified city extents above to compute total light emitted by each city in both years\n","\n","- As a sanity check, make scatterplots of population against night lights (use log scale for better results). Hint: the two should be positevely correlated!"]},{"cell_type":"code","execution_count":null,"id":"031ca5ab-867d-428d-9276-ecfa0c5c6dbf","metadata":{"id":"031ca5ab-867d-428d-9276-ecfa0c5c6dbf"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"3a0a5ed5-3843-4526-b40f-427c561228ba","metadata":{"id":"3a0a5ed5-3843-4526-b40f-427c561228ba"},"source":["Please make scatterplots of population against night lights (in log scale) below. "]},{"cell_type":"code","execution_count":null,"id":"bfb6714d-0434-4e7f-8c50-74b96cb3b453","metadata":{"id":"bfb6714d-0434-4e7f-8c50-74b96cb3b453"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"eadfada7-9e28-4024-9c62-fe8acf5b23ea","metadata":{"id":"eadfada7-9e28-4024-9c62-fe8acf5b23ea"},"source":["# 3. Identify Tanzania's primary city\n","\n","To make things simple, we call the primary city the one with the largest population. \n","\n","- Find the primary city.\n","\n","- Plot the primary city with one color\n","- Plot all cities except the primary with a different color\n","\n","- Add the outline of Tanzania to the plot (this is available in the file `tza_poly/tza_polygon.shp` we have provided for you)"]},{"cell_type":"code","execution_count":null,"id":"d99bf7ca-9e73-4efc-8c86-d3cd84134b57","metadata":{"id":"d99bf7ca-9e73-4efc-8c86-d3cd84134b57"},"outputs":[],"source":[]},{"cell_type":"code","execution_count":null,"id":"20bff3ab-5fe3-4e95-90d4-6d3cf47a14a1","metadata":{"id":"20bff3ab-5fe3-4e95-90d4-6d3cf47a14a1"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"725a79e2-544c-4ec0-abaa-6d9fe88db0bb","metadata":{"id":"725a79e2-544c-4ec0-abaa-6d9fe88db0bb"},"source":["# 4. Create a routable road network\n","\n","As part of the identification stratey in Storeygard (2016), we need to compute the distance between each city and the primary city following the road network. In this section we will prepare the network, in section 5 we will solve it with the Djkstra's algorithm for the least costly path. \n","\n","The network should have the following elements: i) city points (as nodes), ii) roads (as edges), iii) connectors from city to roads (as edges)\n","\n","- Using the file `roads.shp` we have provided for you as well as the point locations of identified cities in Section 1, create a network dataset that you can use to find shortest paths\n","\n","- While this is somewhat involved, you have seen this exact procedure in the lecture on Donaldson and Hornbeck's paper on the American railroad. The steps we show there are what is needed here: you just need to adapt the code to fit the inputs for this application!\n","\n","- For the costs, you can assume the following parameters (these relative costs are made up):\n","    - `cost_per_km_paved` = 1\n","    - `cost_per_km_unpaved` = 4\n","    \n","- For the connector pieces connecting city points to the road network, you can assume that they are paved"]},{"cell_type":"markdown","id":"05e97817-3947-48da-a626-c1785053fdd8","metadata":{"id":"05e97817-3947-48da-a626-c1785053fdd8"},"source":["## 4.1 Create Road connectors to the cities points\n","\n","In this subsection, you will generate the connectors from the cities to the closest road segment"]},{"cell_type":"code","execution_count":null,"id":"4f53f123-ca6a-4d3a-9302-0d5f7042de4f","metadata":{"id":"4f53f123-ca6a-4d3a-9302-0d5f7042de4f"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"a00e51a0-8d77-49c0-b3ac-5433fa9f480c","metadata":{"id":"a00e51a0-8d77-49c0-b3ac-5433fa9f480c"},"source":["## 4.2 Deal with connectivity issues\n","\n","As discussed in class, drawing the connectors is not sufficient to ensure full connectivity due to numerical precision. \n","\n","- Rescale the connectors by a small factor. Hint: There is one city connector that needs an extra large rescaling!\n"]},{"cell_type":"code","execution_count":null,"id":"c5f66d74-a232-4ea7-a6d8-e659a8662151","metadata":{"id":"c5f66d74-a232-4ea7-a6d8-e659a8662151"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"9080e63e-9645-44cb-95fc-8d62d00d2e1b","metadata":{"id":"9080e63e-9645-44cb-95fc-8d62d00d2e1b"},"source":["Remember that after you applied the rescaling you also need to snap the cities' location to the rescaled connectors. This is achieved by the function below"]},{"cell_type":"code","execution_count":null,"id":"1cb45f2b-03f7-4fac-8df3-afbbd8074a2e","metadata":{"id":"1cb45f2b-03f7-4fac-8df3-afbbd8074a2e"},"outputs":[],"source":["def snap_points_to_modified_connectors(points_gdf, connectors_gdf):\n","    \n","    out = points_gdf.copy()\n","    \n","    for index, row in points_gdf.iterrows():\n","        \n","        tmp_gdf = connectors_gdf.copy()\n","        tmp_gdf['distance'] = tmp_gdf.distance(row['geometry'])\n","        closest_geom = list(tmp_gdf.sort_values('distance')['geometry'])[0]\n","        closest_pts = [Point(coords) for coords in closest_geom.coords]\n","        dists = [row['geometry'].distance(x) for x in closest_pts]\n","        minimum = min(dists)\n","        minimum_index = [i for i, j in enumerate(dists) if j == minimum][0]      \n","        out.loc[index, 'geometry'] = closest_pts[minimum_index]\n","        \n","    return out"]},{"cell_type":"markdown","id":"9fff6684-309f-44eb-9894-3283d9a1b744","metadata":{"id":"9fff6684-309f-44eb-9894-3283d9a1b744"},"source":["Use the function above, to snap the identified city to the rescale connectors"]},{"cell_type":"code","execution_count":null,"id":"6be8b211-534e-4b82-8042-337238570288","metadata":{"id":"6be8b211-534e-4b82-8042-337238570288"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"0d14047e-7636-40bb-9d9e-bc616c37ef68","metadata":{"id":"0d14047e-7636-40bb-9d9e-bc616c37ef68"},"source":["## 4.3 Join roads and connectors\n","\n","Now you can join the roads and the rescaled connectors to finalize the construction of your network edges"]},{"cell_type":"code","execution_count":null,"id":"f092cd5a-8ee9-4aa6-a290-2e525d358872","metadata":{"id":"f092cd5a-8ee9-4aa6-a290-2e525d358872"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"2ecf4853-75ab-4983-a3cb-999afb60042f","metadata":{"id":"2ecf4853-75ab-4983-a3cb-999afb60042f"},"source":["## 4.4 Examine connectivity\n","\n","Now we can examine the connectivity of the network you created. \n","\n","Remember that to build the network you need to define the following function"]},{"cell_type":"code","execution_count":null,"id":"155dbdff-461f-44a4-bc80-608541b395c4","metadata":{"id":"155dbdff-461f-44a4-bc80-608541b395c4"},"outputs":[],"source":["# taken from here:\n","# https://www.reddit.com/r/gis/comments/b1ui7h/geopandas_how_to_make_a_graph_out_of_a/\n","\n","def gdf_to_nx(gdf_network):\n","    # generate graph from GeoDataFrame of LineStrings\n","    net = nx.Graph()\n","    net.graph['crs'] = gdf_network.crs\n","    fields = list(gdf_network.columns)\n","\n","    for index, row in gdf_network.iterrows():\n","        first = row.geometry.coords[0]\n","        last = row.geometry.coords[-1]\n","\n","        data = [row[f] for f in fields]\n","        attributes = dict(zip(fields, data))\n","        net.add_edge(first, last, **attributes)\n","\n","    return net"]},{"cell_type":"markdown","id":"d407179d-09bd-494f-9416-fcf3afe3a6f9","metadata":{"id":"d407179d-09bd-494f-9416-fcf3afe3a6f9"},"source":["Build the network!\n","\n","Also generate a geodataframe with the number of connected component for each city. [again check carefully the network notebook for Donalson and Hornbeck's replication]"]},{"cell_type":"code","execution_count":null,"id":"3176af00-e1f9-4598-b186-403544276158","metadata":{"tags":[],"id":"3176af00-e1f9-4598-b186-403544276158"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"2c9cb0e6-687b-4490-bda2-9b93963588c0","metadata":{"id":"2c9cb0e6-687b-4490-bda2-9b93963588c0"},"source":["Now plot it! \n","\n","If you see an unconnected city (number of components=2), then go back to Section 4.2 and inspect the data better!"]},{"cell_type":"code","execution_count":null,"id":"cf49c492-9e7e-4c37-af53-5b3d7bf7fe15","metadata":{"id":"cf49c492-9e7e-4c37-af53-5b3d7bf7fe15"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"48a8cbce-4cc4-4491-b79c-115d7d1d379c","metadata":{"id":"48a8cbce-4cc4-4491-b79c-115d7d1d379c"},"source":["## 4.5 Add costs\n","\n","Now you are ready to assign the costs to paved and unpaved roads."]},{"cell_type":"code","execution_count":null,"id":"35d65eca-eb42-42ac-9616-e9defc847b84","metadata":{"id":"35d65eca-eb42-42ac-9616-e9defc847b84"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"34913f41-7cbc-4d74-86ba-1b4f13a8faa2","metadata":{"id":"34913f41-7cbc-4d74-86ba-1b4f13a8faa2"},"source":["# 5. Find optimal paths between cities\n","\n","Now that your network is routable, fully connected, and costs are specified, you can solve it!\n","\n","- Using the roads and cities network created in step 4, solve the network to find the shortest paths\n","- Once again, you should follow the code we have provided during the lecture on Donaldson and Hornbeck's paper on the American railroad"]},{"cell_type":"code","execution_count":null,"id":"2ba3fc5a-d2b0-4c77-8592-0435e14cbda2","metadata":{"id":"2ba3fc5a-d2b0-4c77-8592-0435e14cbda2"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"06809d50-220a-4382-a332-404ec179b6c9","metadata":{"id":"06809d50-220a-4382-a332-404ec179b6c9"},"source":["- Now iterate over the list of centroids to grab the elements from the dictionary and insert them into a list of lists into the right order. \n","- Finally, convert the  list of lists into a `pandas` dataframe. "]},{"cell_type":"code","execution_count":null,"id":"21759d3e-e4ef-41a5-9e4e-3910595ad9c6","metadata":{"id":"21759d3e-e4ef-41a5-9e4e-3910595ad9c6"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"7aee500d-dfb8-4923-84f0-14fd78980386","metadata":{"id":"7aee500d-dfb8-4923-84f0-14fd78980386"},"source":["Notice that for this application we would just need the distance of each city to the primary city!\n","\n","- Create a dataframe having 2 columns: i) city names and ii) distance to primary city."]},{"cell_type":"code","execution_count":null,"id":"86e78a09-eed1-494f-bf5c-a664f047778b","metadata":{"id":"86e78a09-eed1-494f-bf5c-a664f047778b"},"outputs":[],"source":[]},{"cell_type":"markdown","id":"3be22070-5472-4462-a185-05b6d35cc334","metadata":{"id":"3be22070-5472-4462-a185-05b6d35cc334"},"source":["# 6. Build the final dataset\n","\n","- The point of this data construction exercise is to build a dataset you can use to run Storeygard's regression.\n","\n","- You should therefore build this dataset as a final exercise\n","\n","- The dataset should have\n","    - city names\n","    - city populations\n","    - a column for the year\n","    - lights in both years\n","    - oil prices in both years (we have provided the oil price data in `Europe_Brent_Spot_Price_FOB.csv`)\n","    - distance to the primary city"]},{"cell_type":"code","execution_count":null,"id":"d020b8a5-a404-4590-8f90-36d24bdc166f","metadata":{"id":"d020b8a5-a404-4590-8f90-36d24bdc166f"},"outputs":[],"source":[]}],"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","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.10.4"},"colab":{"provenance":[{"file_id":"1CocWT3ayoySO8krpB0yZQ8FwQwkFzWTj","timestamp":1685891676396},{"file_id":"1jqOgopJ2HQIzCJodAci8B1hrpgrHxHso","timestamp":1684892159946},{"file_id":"17NBn_jJ91p9JV_O0i-6HIjkRHPmpFo3l","timestamp":1684864889423}]}},"nbformat":4,"nbformat_minor":5}