660 lines
695 KiB
Plaintext
660 lines
695 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"id": "initial_id",
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": true,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:10.475296Z",
|
||
|
|
"start_time": "2024-06-10T13:21:10.462294Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"import json\n",
|
||
|
|
"from pathlib import Path\n",
|
||
|
|
"import geopandas as gpd\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"from shapely.geometry import Polygon, Point, LineString, shape\n",
|
||
|
|
"import networkx as nx"
|
||
|
|
],
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 22
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:11.059145Z",
|
||
|
|
"start_time": "2024-06-10T13:21:11.050146Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"buildings_file = \"./input_files/buildings.geojson\"\n",
|
||
|
|
"roads_file = \"./input_files/roads.json\""
|
||
|
|
],
|
||
|
|
"id": "54ab8fc586f52a30",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 23
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:11.589663Z",
|
||
|
|
"start_time": "2024-06-10T13:21:11.571663Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"with open(buildings_file, 'r') as file:\n",
|
||
|
|
" city = json.load(file)"
|
||
|
|
],
|
||
|
|
"id": "1e89883e8ef6043e",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 24
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:12.105043Z",
|
||
|
|
"start_time": "2024-06-10T13:21:12.088015Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"centroids = []\n",
|
||
|
|
"building_ids = [] # List to store building IDs\n",
|
||
|
|
"buildings = city['features']\n",
|
||
|
|
"for building in buildings:\n",
|
||
|
|
" coordinates = building['geometry']['coordinates'][0]\n",
|
||
|
|
" building_polygon = Polygon(coordinates)\n",
|
||
|
|
" centroid = building_polygon.centroid\n",
|
||
|
|
" centroids.append(centroid)\n",
|
||
|
|
" building_ids.append(building['id']) # Extract building ID"
|
||
|
|
],
|
||
|
|
"id": "43231d3a9c2831c",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 25
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:14.131575Z",
|
||
|
|
"start_time": "2024-06-10T13:21:14.118038Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"import json\n",
|
||
|
|
"from pathlib import Path\n",
|
||
|
|
"from shapely.geometry import Polygon, Point, shape\n",
|
||
|
|
"\n",
|
||
|
|
"def process_geojson(x, y, diff):\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" Processes a GeoJSON file to find roads that have at least one node within a specified polygon.\n",
|
||
|
|
" \n",
|
||
|
|
" Parameters:\n",
|
||
|
|
" x (float): The x-coordinate of the center of the selection box.\n",
|
||
|
|
" y (float): The y-coordinate of the center of the selection box.\n",
|
||
|
|
" diff (float): The half-width of the selection box.\n",
|
||
|
|
" \n",
|
||
|
|
" Returns:\n",
|
||
|
|
" str: The file path of the output GeoJSON file containing the selected roads.\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" diff += 2*diff\n",
|
||
|
|
" # Define the selection box (polygon)\n",
|
||
|
|
" selection_box = Polygon([\n",
|
||
|
|
" [x + diff, y - diff],\n",
|
||
|
|
" [x - diff, y - diff],\n",
|
||
|
|
" [x - diff, y + diff],\n",
|
||
|
|
" [x + diff, y + diff]\n",
|
||
|
|
" ])\n",
|
||
|
|
"\n",
|
||
|
|
" # Define input and output file paths\n",
|
||
|
|
" geojson_file = Path(\"./input_files/roads.json\").resolve()\n",
|
||
|
|
" output_file = Path('./input_files/output_roads.geojson').resolve()\n",
|
||
|
|
" \n",
|
||
|
|
" # Initialize a list to store the roads in the region\n",
|
||
|
|
" roads_in_region = []\n",
|
||
|
|
"\n",
|
||
|
|
" # Read the GeoJSON file\n",
|
||
|
|
" with open(geojson_file, 'r') as file:\n",
|
||
|
|
" roads = json.load(file)\n",
|
||
|
|
" line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']\n",
|
||
|
|
"\n",
|
||
|
|
" # Check each road feature\n",
|
||
|
|
" for feature in line_features:\n",
|
||
|
|
" road_shape = shape(feature['geometry'])\n",
|
||
|
|
" # Check if any node of the road is inside the selection box\n",
|
||
|
|
" if any(selection_box.contains(Point(coord)) for coord in road_shape.coords):\n",
|
||
|
|
" roads_in_region.append(feature)\n",
|
||
|
|
"\n",
|
||
|
|
" # Create a new GeoJSON structure with the selected roads\n",
|
||
|
|
" output_geojson = {\n",
|
||
|
|
" \"type\": \"FeatureCollection\",\n",
|
||
|
|
" \"features\": roads_in_region\n",
|
||
|
|
" }\n",
|
||
|
|
"\n",
|
||
|
|
" # Write the selected roads to the output file\n",
|
||
|
|
" with open(output_file, 'w') as outfile:\n",
|
||
|
|
" json.dump(output_geojson, outfile)\n",
|
||
|
|
"\n",
|
||
|
|
" return output_file"
|
||
|
|
],
|
||
|
|
"id": "26792dbde788f8b6",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 26
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:21.767454Z",
|
||
|
|
"start_time": "2024-06-10T13:21:14.891727Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": "road_geojson = process_geojson(-73.61038745537758, 45.484399882086215, 0.001)",
|
||
|
|
"id": "ec7c16b59965fdf",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 27
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:23.297352Z",
|
||
|
|
"start_time": "2024-06-10T13:21:23.283323Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"with open(road_geojson, 'r') as file:\n",
|
||
|
|
" roads = json.load(file)"
|
||
|
|
],
|
||
|
|
"id": "746cf9f190d6390a",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 28
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:23.801263Z",
|
||
|
|
"start_time": "2024-06-10T13:21:23.789266Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"line_features = [feature for feature in roads['features'] if feature['geometry']['type'] == 'LineString']\n",
|
||
|
|
"\n",
|
||
|
|
"# Create a list of LineString objects and their properties\n",
|
||
|
|
"lines = []\n",
|
||
|
|
"for feature in line_features:\n",
|
||
|
|
" # Create a LineString from coordinates\n",
|
||
|
|
" linestring = LineString(feature['geometry']['coordinates'])\n",
|
||
|
|
" lines.append(linestring)"
|
||
|
|
],
|
||
|
|
"id": "9dfe16b2d4d98333",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 29
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:24.368487Z",
|
||
|
|
"start_time": "2024-06-10T13:21:24.349490Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"cleaned_lines = []\n",
|
||
|
|
"for line in lines:\n",
|
||
|
|
" coords = list(line.coords)\n",
|
||
|
|
" cleaned_line = LineString([coords[0], coords[-1]])\n",
|
||
|
|
" cleaned_lines.append(cleaned_line)"
|
||
|
|
],
|
||
|
|
"id": "c15a91fd07f8e25e",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 30
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:25.169119Z",
|
||
|
|
"start_time": "2024-06-10T13:21:24.890980Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"fig, ax = plt.subplots()\n",
|
||
|
|
" \n",
|
||
|
|
"for line in cleaned_lines:\n",
|
||
|
|
" # Extract the x and y coordinates of the LineString\n",
|
||
|
|
" x, y = line.xy\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot the LineString\n",
|
||
|
|
" ax.plot(x, y, color='blue')\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot the nodes (points)\n",
|
||
|
|
" for coord in line.coords:\n",
|
||
|
|
" point = Point(coord)\n",
|
||
|
|
" ax.plot(point.x, point.y, 'ro') # 'ro' for red circles\n",
|
||
|
|
"\n",
|
||
|
|
"ax.set_title('LineString Objects and Their Nodes')\n",
|
||
|
|
"ax.set_xlabel('Longitude')\n",
|
||
|
|
"ax.set_ylabel('Latitude')\n",
|
||
|
|
"plt.show()"
|
||
|
|
],
|
||
|
|
"id": "5d8fbd8f6bbae402",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 640x480 with 1 Axes>"
|
||
|
|
],
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAHHCAYAAABEEKc/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADEnklEQVR4nOydd1hT5xfHv2EFRIYDAQUB967VCrVuUakTnHXU1aq4qrbWVn9q1S6t2qp1gXvv2bpH1da9N+6FCjgZyk7O749jEgIBI5LchLyf58kDvHkTTu7Nvfd7z3uGjIgIAoFAIBAIBAI1VlIbIBAIBAKBQGBqCIEkEAgEAoFAkAkhkAQCgUAgEAgyIQSSQCAQCAQCQSaEQBIIBAKBQCDIhBBIAoFAIBAIBJkQAkkgEAgEAoEgE0IgCQQCgUAgEGRCCCSBQCAQCASCTAiBJLBI7t27B5lMhiVLlkhtynvh6+uLXr16SW2GXqi2+dSpU986d/z48ZDJZEawyvzI622jer9nz57l2Xv26tULvr6+efZ+xmDJkiWQyWS4d++e1KYITAQhkAT5DtWJ7vTp01KbAqVSiWXLliEgIACFCxeGk5MTypUrhx49euD48ePqeVevXsX48ePN7uT8+vVr/PTTT6hWrRoKFCgAFxcX1KtXD8uWLYMpdzHasWMHxo8fL7UZBuHgwYOQyWR6PUwZX19fyGQyfPXVV1meU33GDRs2SGCZwFKwkdoAgUAKfHx8kJSUBFtbW4P+nyFDhmD27NkIDg5Gt27dYGNjg+vXr2Pnzp0oVaoUPv74YwAskCZMmICGDRu+05339evXYWUlzX1OTEwMAgMDERERgc6dO2Pw4MFITk7Gxo0b0bNnT+zYsQMrV66EtbX1O7/3mDFjMHLkSANYzezYsQOzZ8/OlyKpYsWKWL58udbYqFGjULBgQYwePdooNsyfPx9KpTLP3mvUqFEoXrx4nryfQKAvQiAJLBKZTAZ7e3uD/o+YmBjMmTMHffv2xbx587Semz59Op4+fZqr9yUiJCcnw8HBAXK5PC9MzRU9e/ZEREQENm/ejDZt2qjHhwwZghEjRmDq1Kn48MMP8f3337/ze9vY2MDGRpyecoO7uzs+//xzrbFJkyahaNGiWcYNhT43Hunp6VAqlbCzs8t2TuXKlXH9+nVMmjQJf/75Z16aKBC8FbHEJrBIdMUg9erVCwULFsSjR48QEhKCggULws3NDd9++y0UCoXW65VKJaZPn47KlSvD3t4e7u7uCA0NxcuXL9Vz7t69CyJCnTp1svx/mUyGYsWKAeAlwY4dOwIAGjVqpF7+OHjwIABeamjVqhV2796Njz76CA4ODggPD1c/lzEGSbW8eOTIEXzzzTdwc3ODo6Mj2rZtm0WQKZVKjB8/HsWLF0eBAgXQqFEjXL16Va+4puPHj2P37t3o1auXljhSMXHiRJQtWxa//fYbkpKSsjw/bdo0+Pj4wMHBAQ0aNMDly5e1ns8uzmbFihWoWbMmHBwcULhwYXTu3BmRkZFZ5p04cQItWrRAoUKF4OjoiGrVqmHGjBkAeD/Pnj0bAHQuN61ZswY1a9aEk5MTnJ2dUbVqVfVrc2Lq1Kn45JNPUKRIETg4OKBmzZo6l4BkMhkGDx6MLVu2oEqVKpDL5ahcuTJ27dqVZe7hw4dRq1Yt2Nvbo3Tp0ur9bghiY2PRq1cvuLq6wsXFBb1790ZiYmKWefrsg8wxSBnjz6ZPn47SpUtDLpfj6tWrOdrk6+uLHj16YP78+Xj8+PFbP8O5c+fQvHlzODs7o2DBgggMDNRaylZx5coVNG7cGA4ODvDy8sLPP/+crcdr586dqFevHhwdHeHk5ISWLVviypUrWnOio6PRu3dveHl5QS6Xw9PTE8HBwWa3ZC7QRtyiCQQZUCgUCAoKQkBAAKZOnYp9+/bh999/R+nSpTFgwAD1vNDQUCxZsgS9e/fGkCFDcPfuXcyaNQvnzp3DkSNHYGtrCx8fHwDA+vXr0bFjRxQoUEDn/6xfvz6GDBmCP//8E//73/9QsWJFAFD/BHgprUuXLggNDUXfvn1Rvnz5HD/HV199hUKFCmHcuHG4d+8epk+fjsGDB2Pt2rXqOaNGjcLkyZPRunVrBAUF4cKFCwgKCkJycvJbt9Pff/8NAOjRo4fO521sbNC1a1dMmDABR44cQZMmTdTPLVu2DAkJCRg0aBCSk5MxY8YMNG7cGJcuXYK7u3u2//OXX37B2LFj0alTJ/Tp0wdPnz7FzJkzUb9+fZw7dw6urq4AgL1796JVq1bw9PTE0KFD4eHhgYiICGzbtg1Dhw5FaGgoHj9+jL1792ZZitq7dy+6dOmCwMBA/PbbbwCAiIgIHDlyBEOHDs1xm8yYMQNt2rRBt27dkJqaijVr1qBjx47Ytm0bWrZsqTX38OHD2LRpEwYOHAgnJyf8+eefaN++PR48eIAiRYoAAC5duoRmzZrBzc0N48ePR3p6OsaNG5fjNnofOnXqBD8/P0ycOBFnz57FggULUKxYMfV2APTfB9mxePFiJCcno1+/fpDL5ShcuPBb7Ro9ejSWLVv2Vi/SlStXUK9ePTg7O+O7776Dra0twsPD0bBhQxw6dAgBAQEAWMw0atQI6enpGDlyJBwdHTFv3jw4ODhkec/ly5ejZ8+eCAoKwm+//YbExETMnTsXdevWxblz59QisH379rhy5Qq++uor+Pr64smTJ9i7dy8ePHhgdsHqggyQQJDPWLx4MQGgU6dOZTvn7t27BIAWL16sHuvZsycBoB9//FFr7ocffkg1a9ZU//3ff/8RAFq5cqXWvF27dmUZ79GjBwGgQoUKUdu2bWnq1KkUERGRxZ7169cTADpw4ECW53x8fAgA7dq1S+dzPXv2zPLZmzRpQkqlUj3+9ddfk7W1NcXGxhIRUXR0NNnY2FBISIjW+40fP54AaL2nLkJCQggAvXz5Mts5mzZtIgD0559/EpFmmzs4ONDDhw/V806cOEEA6Ouvv1aPjRs3jjKenu7du0fW1tb0yy+/aP2PS5cukY2NjXo8PT2d/Pz8yMfHJ4ttGbfHoEGDSNfpb+jQoeTs7Ezp6ek5fn5dJCYmav2dmppKVapUocaNG2uNAyA7Ozu6deuWeuzChQsEgGbOnKkeCwkJIXt7e7p//7567OrVq2Rtba3T9pyoXLkyNWjQQOdzqm39xRdfaI23bduWihQpov5b331AxMeSj4+P+m/Vvnd2dqYnT57oZbOPjw+1bNmSiIh69+5N9vb29PjxYyIiOnDgAAGg9evXq+eHhISQnZ0d3b59Wz32+PFjcnJyovr166vHhg0bRgDoxIkT6rEnT56Qi4sLAaC7d+8SEVFCQgK5urpS3759teyKjo4mFxcX9fjLly8JAE2ZMkWvzyUwH8QSm0CQif79+2v9Xa9ePdy5c0f99/r16+Hi4oKmTZvi2bNn6kfNmjVRsGBBHDhwQD138eLFmDVrFvz8/LB582Z8++23qFixIgIDA/Ho0SO9bfLz80NQUJDe8/v166e1bFSvXj0oFArcv38fALB//36kp6dj4MCBWq/TlTGki4SEBACAk5NTtnNUz8XHx2uNh4SEoESJEuq//f39ERAQgB07dmT7Xps2bYJSqUSnTp20trmHhwfKli2r3ubnzp3D3bt3MWzYsCzeDH2ytlxdXfH69Wvs3bv3rXMzk9ED8fLlS8TFxaFevXo4e/ZslrlNmjRB6dKl1X9Xq1YNzs7O6u+ZQqHA7t27ERISgpIlS6rnVaxY8Z2+B++Cru/98+fP1ftP332QE+3bt4ebm9s72zZmzBikp6dj0qRJOp9XKBTYs2cPQkJCUKpUKfW4p6cnunbtisOHD6s/x44dO/Dxxx/D399fPc/NzQ3dunXTes+9e/ciNjYWXbp00fq81tbWCAgIUH9eBwcH2NnZ4eDBg1pL7ALzRwgkCUhJSUH16tUhk8lw/vz5HOc2bNgwS2pu5hOZiufPn8PLywsymQyxsbFaz61cuRIffPABChQoAE9PT3zxxRd4/vz5O9n9yy+/4JNPPkGBAgXe6ko3V+zt7bOcwAsVKqR
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 31
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:26.302823Z",
|
||
|
|
"start_time": "2024-06-10T13:21:26.267819Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"closest_roads = []\n",
|
||
|
|
"unique_roads_set = set()\n",
|
||
|
|
"\n",
|
||
|
|
"# Loop through each centroid\n",
|
||
|
|
"for centroid in centroids:\n",
|
||
|
|
" min_distance = float('inf') # Start with a large number to ensure any real distance is smaller\n",
|
||
|
|
" closest_road = None\n",
|
||
|
|
" \n",
|
||
|
|
" # Loop through each road and calculate the distance to the current centroid\n",
|
||
|
|
" for line in cleaned_lines:\n",
|
||
|
|
" distance = line.distance(centroid)\n",
|
||
|
|
" # Check if the current road is closer than the ones previously checked\n",
|
||
|
|
" if distance < min_distance:\n",
|
||
|
|
" min_distance = distance\n",
|
||
|
|
" closest_road = line\n",
|
||
|
|
"\n",
|
||
|
|
" # Add the closest road to the list if it's not already added\n",
|
||
|
|
" if closest_road and closest_road.wkt not in unique_roads_set:\n",
|
||
|
|
" unique_roads_set.add(closest_road.wkt)\n",
|
||
|
|
" closest_roads.append(closest_road)"
|
||
|
|
],
|
||
|
|
"id": "a2ca0b96ca869007",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 32
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:27.461902Z",
|
||
|
|
"start_time": "2024-06-10T13:21:27.246748Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"fig, ax = plt.subplots()\n",
|
||
|
|
" \n",
|
||
|
|
"for line in closest_roads:\n",
|
||
|
|
" # Extract the x and y coordinates of the LineString\n",
|
||
|
|
" x, y = line.xy\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot the LineString\n",
|
||
|
|
" ax.plot(x, y, color='blue')\n",
|
||
|
|
" \n",
|
||
|
|
" # Plot the nodes (points)\n",
|
||
|
|
" for coord in line.coords:\n",
|
||
|
|
" point = Point(coord)\n",
|
||
|
|
" ax.plot(point.x, point.y, 'ro') # 'ro' for red circles\n",
|
||
|
|
"\n",
|
||
|
|
"ax.set_title('LineString Objects and Their Nodes')\n",
|
||
|
|
"ax.set_xlabel('Longitude')\n",
|
||
|
|
"ax.set_ylabel('Latitude')"
|
||
|
|
],
|
||
|
|
"id": "8a358c26192dc7a0",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"Text(0, 0.5, 'Latitude')"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 33,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 640x480 with 1 Axes>"
|
||
|
|
],
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHHCAYAAACfqw0dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACVY0lEQVR4nOzdeXhM5xcH8O9kkkwWkiArQkIRewhJQ+whbRWpttS+1VJalNpqCYrEVltR2tqXEJTWHopagpaIJRFK1JrYE7Ukkjm/P84vw8giicTNJOfzPPOkc+edmTPXNPfkXc6rIiKCEEIIIYTIFiOlAxBCCCGEMESSRAkhhBBC5IAkUUIIIYQQOSBJlBBCCCFEDkgSJYQQQgiRA5JECSGEEELkgCRRQgghhBA5IEmUEEIIIUQOSBIlhBBCCJEDkkQJkYErV65ApVJh2bJlSofyRlxcXNC9e3elw8iS1HM+Y8aM17YdP348VCrVW4jK8OT2uUl9vbt37+baa3bv3h0uLi659npvw7Jly6BSqXDlyhWlQxH5hCRRolBK/WX4999/Kx0KtFotVqxYAS8vLxQvXhxFixZFxYoV0bVrVxw9elTXLjIyEuPHjze4X+CPHz/Gd999hxo1asDCwgLW1tZo0KABVqxYgfy869T27dsxfvx4pcPIE/v374dKpcrSLT9zcXGBSqXCV199leax1M+4YcMGBSIThYWx0gEIkV+VLVsWT58+hYmJSZ6+z8CBAzF//ny0adMGnTp1grGxMaKjo7Fjxw6UK1cO7777LgBOoiZMmIDGjRtn6y/46OhoGBkp8/dSXFwcmjVrhqioKHz22Wf48ssv8ezZM2zcuBHdunXD9u3bsXr1aqjV6my/9pgxYzBy5Mg8iJpt374d8+fPL5CJVOXKlbFy5Uq9Y6NGjUKRIkUwevTotxLDTz/9BK1Wm2uvNWrUKJQsWTJXXk+IrJIkSogMqFQqmJmZ5el7xMXFYcGCBejduzcWL16s99js2bNx586dHL0uEeHZs2cwNzeHRqPJjVBzpFu3boiKisKvv/6K1q1b644PHDgQw4YNw4wZM1CrVi2MGDEi269tbGwMY2P5FZYTDg4O6Ny5s96xoKAg2NrapjmeV7Lyx0lycjK0Wi1MTU0zbFO1alVER0cjKCgIc+fOzc0QhXgtGc4TIgPpzYnq3r07ihQpghs3bsDf3x9FihSBnZ0dvvnmG6SkpOg9X6vVYvbs2ahatSrMzMzg4OCAvn374sGDB7o2MTExICLUr18/zfurVCrY29sD4OHHTz/9FADQpEkT3VDL/v37AfCwxocffohdu3ahTp06MDc3x6JFi3SPvTwnKnUo8/DhwxgyZAjs7OxgaWmJjz76KE3SptVqMX78eJQsWRIWFhZo0qQJIiMjszTP6ujRo9i1axe6d++ul0ClCgwMRIUKFTB16lQ8ffo0zeOzZs1C2bJlYW5ujkaNGuHs2bN6j2c072fVqlXw8PCAubk5ihcvjs8++wzXrl1L0+7YsWP44IMPUKxYMVhaWqJGjRqYM2cOAP53nj9/PgCkO7QVHBwMDw8PFC1aFFZWVqhevbruuZmZMWMG6tWrhxIlSsDc3BweHh7pDjepVCp8+eWX2Lx5M6pVqwaNRoOqVati586dadoeOnQIdevWhZmZGcqXL6/7d88LDx8+RPfu3WFjYwNra2v06NEDT548SdMuK/8Gr86Jenk+3OzZs1G+fHloNBpERkZmGpOLiwu6du2Kn376CTdv3nztZwgPD8f7778PKysrFClSBM2aNdMbNk917tw5NG3aFObm5ihdujQmTZqUYc/Zjh070KBBA1haWqJo0aJo2bIlzp07p9cmNjYWPXr0QOnSpaHRaODk5IQ2bdoY3PC80Cd/xgmRTSkpKfDz84OXlxdmzJiBPXv2YObMmShfvjy++OILXbu+ffti2bJl6NGjBwYOHIiYmBj88MMPCA8Px+HDh2FiYoKyZcsCAEJCQvDpp5/CwsIi3fds2LAhBg4ciLlz5+Lbb79F5cqVAUD3E+Bhuw4dOqBv377o3bs3KlWqlOnn+Oqrr1CsWDEEBATgypUrmD17Nr788kusW7dO12bUqFGYNm0aWrVqBT8/P0RERMDPzw/Pnj177Xn6/fffAQBdu3ZN93FjY2N07NgREyZMwOHDh+Hr66t7bMWKFXj06BEGDBiAZ8+eYc6cOWjatCnOnDkDBweHDN9z8uTJGDt2LNq1a4fPP/8cd+7cwbx589CwYUOEh4fDxsYGABAaGooPP/wQTk5OGDRoEBwdHREVFYWtW7di0KBB6Nu3L27evInQ0NA0w16hoaHo0KEDmjVrhqlTpwIAoqKicPjwYQwaNCjTczJnzhy0bt0anTp1QlJSEoKDg/Hpp59i69ataNmypV7bQ4cOYdOmTejfvz+KFi2KuXPn4uOPP8bVq1dRokQJAMCZM2fQokUL2NnZYfz48UhOTkZAQECm5+hNtGvXDq6urggMDMTJkyfx888/w97eXncegKz/G2Rk6dKlePbsGfr06QONRoPixYu/Nq7Ro0djxYoVr+2NOnfuHBo0aAArKysMHz4cJiYmWLRoERo3bowDBw7Ay8sLACc8TZo0QXJyMkaOHAlLS0ssXrwY5ubmaV5z5cqV6NatG/z8/DB16lQ8efIECxcuhI+PD8LDw3WJ4scff4xz587hq6++gouLC27fvo3Q0FBcvXrV4CbYi5eQEIXQ0qVLCQD99ddfGbaJiYkhALR06VLdsW7duhEAmjhxol7bWrVqkYeHh+7+wYMHCQCtXr1ar93OnTvTHO/atSsBoGLFitFHH31EM2bMoKioqDTxhISEEADat29fmsfKli1LAGjnzp3pPtatW7c0n93X15e0Wq3u+Ndff01qtZoePnxIRESxsbFkbGxM/v7+eq83fvx4AqD3munx9/cnAPTgwYMM22zatIkA0Ny5c4noxTk3Nzen69ev69odO3aMANDXX3+tOxYQEEAv/wq7cuUKqdVqmjx5st57nDlzhoyNjXXHk5OTydXVlcqWLZsmtpfPx4ABAyi9X5GDBg0iKysrSk5OzvTzp+fJkyd695OSkqhatWrUtGlTveMAyNTUlP755x/dsYiICAJA8+bN0x3z9/cnMzMz+vfff3XHIiMjSa1Wpxt7ZqpWrUqNGjVK97HUc92zZ0+94x999BGVKFFCdz+r/wZE/P9S2bJldfdT/+2trKzo9u3bWYq5bNmy1LJlSyIi6tGjB5mZmdHNmzeJiGjfvn0EgEJCQnTt/f39ydTUlC5duqQ7dvPmTSpatCg1bNhQd2zw4MEEgI4dO6Y7dvv2bbK2tiYAFBMTQ0REjx49IhsbG+rdu7deXLGxsWRtba07/uDBAwJA06dPz9LnEoZDhvOEyIF+/frp3W/QoAEuX76sux8SEgJra2s0b94cd+/e1d08PDxQpEgR7Nu3T9d26dKl+OGHH+Dq6opff/0V33zzDSpXroxmzZrhxo0bWY7J1dUVfn5+WW7fp08fvSGqBg0aICUlBf/++y8AYO/evUhOTkb//v31npfeSqj0PHr0CABQtGjRDNukPpaQkKB33N/fH6VKldLd9/T0hJeXF7Zv357ha23atAlarRbt2rXTO+eOjo6oUKGC7pyHh4cjJiYGgwcPTtMrkpXVaDY2Nnj8+DFCQ0Nf2/ZVL/dkPHjwAPHx8WjQoAFOnjyZpq2vry/Kly+vu1+jRg1YWVnpvmcpKSnYtWsX/P39UaZMGV27ypUrZ+t7kB3pfe/v3bun+/fL6r9BZj7++GPY2dllO7YxY8YgOTkZQUFB6T6ekpKC3bt3w9/fH+XKldMdd3JyQseOHXHo0CHd59i+fTveffddeHp66trZ2dmhU6dOeq8ZGhqKhw8fokOHDnqfV61Ww8vLS/d5zc3NYWpqiv379+sN5wvDJ0lUPpWYmAh3d3eoVCqcOnUq07aNGzdOsyz51V92qe7du4fSpUtDpVLh4cOHeo+tXr0aNWvWhIWFBZycnNCzZ0/cu3cvW3FPnjwZ9erVg4WFxWu77Q2VmZlZml/yxYoV0/vlePHiRcT
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 33
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:30.801378Z",
|
||
|
|
"start_time": "2024-06-10T13:21:30.787380Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"def find_nearest_point_on_line(point, line):\n",
|
||
|
|
" return line.interpolate(line.project(point))\n",
|
||
|
|
"\n",
|
||
|
|
"# List to store the nearest points\n",
|
||
|
|
"nearest_points = []\n",
|
||
|
|
"\n",
|
||
|
|
"# Find the nearest point on each closest road for each centroid\n",
|
||
|
|
"for centroid in centroids:\n",
|
||
|
|
" # Find the closest road for this centroid\n",
|
||
|
|
" min_distance = float('inf')\n",
|
||
|
|
" closest_road = None\n",
|
||
|
|
" for road in closest_roads:\n",
|
||
|
|
" distance = centroid.distance(road)\n",
|
||
|
|
" if distance < min_distance:\n",
|
||
|
|
" min_distance = distance\n",
|
||
|
|
" closest_road = road\n",
|
||
|
|
"\n",
|
||
|
|
" # Find the nearest point on the closest road\n",
|
||
|
|
" if closest_road:\n",
|
||
|
|
" nearest_point = find_nearest_point_on_line(centroid, closest_road)\n",
|
||
|
|
" nearest_points.append(nearest_point)"
|
||
|
|
],
|
||
|
|
"id": "2cc9ba52c2e08c34",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 34
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:31.668532Z",
|
||
|
|
"start_time": "2024-06-10T13:21:31.659532Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": "len(nearest_points)",
|
||
|
|
"id": "ba24ef6761900721",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"57"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 35,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 35
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:32.579043Z",
|
||
|
|
"start_time": "2024-06-10T13:21:32.554017Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"def break_down_roads(closest_roads, nearest_points_list):\n",
|
||
|
|
" new_segments = []\n",
|
||
|
|
" for road in closest_roads:\n",
|
||
|
|
" # Get coordinates of the road\n",
|
||
|
|
" coords = list(road.coords)\n",
|
||
|
|
" # Find all nearest points for this road\n",
|
||
|
|
" points_on_road = [point for point in nearest_points_list if road.distance(point) < 0.000000001]\n",
|
||
|
|
" # Sort nearest points along the road\n",
|
||
|
|
" sorted_points = sorted(points_on_road, key=lambda point: road.project(point))\n",
|
||
|
|
" # Add the start node to the sorted points\n",
|
||
|
|
" sorted_points.insert(0, Point(coords[0]))\n",
|
||
|
|
" # Add the end node to the sorted points\n",
|
||
|
|
" sorted_points.append(Point(coords[-1]))\n",
|
||
|
|
" # Create new segments\n",
|
||
|
|
" for i in range(len(sorted_points) - 1):\n",
|
||
|
|
" segment = LineString([sorted_points[i], sorted_points[i + 1]])\n",
|
||
|
|
" new_segments.append(segment)\n",
|
||
|
|
" return new_segments\n",
|
||
|
|
"\n",
|
||
|
|
"# Create new segments\n",
|
||
|
|
"new_segments = break_down_roads(closest_roads, nearest_points)\n",
|
||
|
|
"\n",
|
||
|
|
"len(new_segments)"
|
||
|
|
],
|
||
|
|
"id": "3766eaeb084545f9",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"61"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 36,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 36
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:33.783009Z",
|
||
|
|
"start_time": "2024-06-10T13:21:33.554797Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"import random\n",
|
||
|
|
"def random_color():\n",
|
||
|
|
" return \"#\"+''.join([random.choice('0123456789ABCDEF') for j in range(6)])\n",
|
||
|
|
"fig, ax = plt.subplots(figsize=(15, 10)) # Increase the figure size\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot the segments with random colors\n",
|
||
|
|
"for i, segment in enumerate(new_segments):\n",
|
||
|
|
" x, y = segment.xy\n",
|
||
|
|
" ax.plot(x, y, color=random_color(), linewidth=3, label=f'Segment {i+1}') # Increase the linewidt\n",
|
||
|
|
"\n",
|
||
|
|
"ax.set_title('Road Segments and Nearest Points', fontsize=16)\n",
|
||
|
|
"ax.set_xlabel('Longitude', fontsize=14)\n",
|
||
|
|
"ax.set_ylabel('Latitude', fontsize=14)\n",
|
||
|
|
"plt.show()"
|
||
|
|
],
|
||
|
|
"id": "ac8694fc95c5831e",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 1500x1000 with 1 Axes>"
|
||
|
|
],
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABPEAAANeCAYAAABzsJstAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdZ3gUdfv28XN30wNJIEBCaAldpSNNUKqEIgLSQYiAoqIoYv/bvW/FXgGx01SKigVBuqKCIF1pgiT00CGQnt3f88KH3MbMYkDYTML3cxw5kJlrdq/ZLOvuuTPXOIwxRgAAAAAAAABsy1nYDQAAAAAAAAA4O0I8AAAAAAAAwOYI8QAAAAAAAACbI8QDAAAAAAAAbI4QDwAAAAAAALA5QjwAAAAAAADA5gjxAAAAAAAAAJsjxAMAAAAAAABsjhAPAAAAAAAAsDlCPAAACkFsbKwcDkeen8DAQFWsWFHdu3fXnDlzCrvFf/Tdd9/J4XCoTZs257zt2rVrNWzYMFWvXl3BwcEKCQlRlSpV1LJlS913331auHDhhW8Yl5Qnn3xSDodDTz75ZIG3OfOcdjgcio6OVmpqqmXd3r17c+vw79100035Xg/9/PxUtmxZXXvttZoyZYqMMf/6fiZNmiSHw6Gbbrrp3zcNAEAhIMQDAKAQtWzZUgkJCUpISFCXLl3k5+enr776St26ddOYMWMKu72L4s0331STJk304YcfKiMjQ23btlXPnj1Vq1Ytbdu2TS+//LIefvjhwm6zWGjTpo0cDoe+++67wm6lyDl48KBefvnlwm6jSDifwNRKtWrVcl8Pe/furaioKC1atCj37263+8I0fAER5gIAfMmvsBsAAOBSdvPNN+c5KiQnJ0f33HOPxo0bp1dffVUDBgxQkyZNCq/BC2zjxo0aPXq0PB6PXn31VY0aNUoulyt3vcfj0Y8//qgff/yxELvEpS44OFgZGRl66aWXdPvtt6ts2bKF3dIloVWrVpo0aVKeZW+99ZZGjhypzz//XJMnT9awYcPO+/Z79uyp5s2bKzw8/F92CgBA4eBIPAAAbMTPz08vvviiwsLCJElff/11IXd0Yc2aNUsej0ctWrTQ6NGj8wR4kuR0OnXNNdfo//7v/wqpQ0CKiYlR7969derUKf33v/8t7HYuabfffrtat24tSZo5c+a/uq3w8HDVrl1b5cuXvxCtAQDgc4R4AADYTFBQkGrUqCHpz1P6/i4nJ0cTJ07UVVddpfDw8Nz6u+66S/v27bO8zVWrVumBBx5Q06ZNFR0drYCAAEVFRalbt25atGjRWfuZMmWKmjRpopCQEJUuXVqdOnXSDz/8cF77dmZ/ypUrd17bp6en6+WXX1bz5s0VERGhoKAg1apVSw888ICOHj1quY0xRh988IGuvPJKhYSEKDIyUp07d9by5cu9zvX76/LMzEw99dRTqlmzpoKCglS5cmU9+OCDysjIkCSdPHlS9913n6pWraqgoCDFxsbqySefVE5Ojtf9WLx4sW644QaVL19eAQEBKleunHr27KkVK1ZY1v/1lL3PPvtMrVq1UlhYmEJDQ9WyZUvNnTvXsv/vv/9ektS2bds888b+erTTmjVr1K9fP1WsWFEBAQEKCwtT1apV1atXL3355Zdn/X38VXZ2tqZNm6ZBgwapdu3aCgsLU3BwsGrVqqW77rpL+/fvt9zur6f8rl+/XjfccIPKlCmjwMBAXX755Xr55Ze9zkNLT0/Xk08+qRo1aigwMFDly5dXQkKCdu/eXeC+vXnmmWfk5+eniRMnKjEx8Zy2zcnJ0Xvvvac2bdqodOnSCgwMVFxcnG6//Xbt2bPHcpvPP/9cN998s+rUqaNSpUopKChIcXFxGjZsmLZt22a5zZlZcpMmTdJvv/2mfv36qXz58nK5XHlObT2ffhYtWqRu3bopKipK/v7+KlWqlGrUqKEbb7xRy5Yty61zOBx66qmnJElPPfVUnufZhZo917hxY0lSUlJSnuVbt27V0KFDVaVKFQUGBqp06dJq376917DP20y8v/57z87O1vPPP68rrrhCwcHBioyM1A033KAtW7bk2ebMKcRn/H2m3197nTVrljp06KDIyEj5+/srMjJSl19+uW655RZt3Ljx/B8YAMClxwAAAJ+rUqWKkWQ+/PBDy/U1atQwksxjjz2WZ3lGRobp0KGDkWSCgoJM586dTb9+/UylSpWMJFOmTBmzZs2afLfXvn1743Q6Td26dU2XLl1Mnz59TKNGjYwkI8m89tprln3cddddRpJxOp3mmmuuMf379zeXX365cTqd5u677zaSTOvWrQu83//5z3+MJFOiRAnz66+/Fng7Y4zZt2+fqVu3rpFkSpcubTp06GB69uyZ+1jGxsaapKSkfNvdfvvtufvQunVr079/f3PFFVcYl8tl7r33Xst9WLp0qZFkWrRoYVq3bm3CwsLM9ddfb6677joTHh5uJJnrrrvOHD161NSqVcuULVvW9OrVy3Ts2NEEBQUZSea2226z3I8z9+l0Ok3Tpk1Nnz59TLNmzYzD4TAul8t88MEH+bY583t6/PHHjcPhMC1btjT9+vUz9evXN5KMw+Ewn3/+eW79li1bTEJCgomKijKSTHx8vElISMj9+eGHH4wxxixatMj4+/sbSaZ+/fqmd+/epmfPnqZp06YmMDDQdO/evcC/nz179hhJJjw83DRv3tz06dPHdOnSxcTExBhJpmzZsmb79u35tmvdurWRZB566CETEBBgLrvsMtO/f3/TunVr43K5jCRz991359suNTXVNG/e3EgyoaGh5rrrrjN9+vQxUVFRJjIy0gwZMsRIMk888USB9+HM771atWrGmP89dwYOHGi5r1ZvpVNSUkybNm1yn+etW7c2vXv3NrVq1TKSTGRkpFm7dm2+7VwulwkJCTFXXnmlueGGG8z1119vqlatmrt/P/30U75tEhISjCRzyy23mMDAQBMbG2v69u1runXrZl566aXz7mfSpEnG4XAYh8NhmjVrZvr162euv/5606hRI+NyufL8PhISEnKfh/Xr18/zPHv33XcL9Lif2Y+EhATL9TfffLORZOrVq5e7bM6cObn/1mrVqmX69+9v2rVrl/ucGTZsWL7b+fDDDy3v58zv/aqrrjIdOnQwISEhplOnTqZXr165r60REREmMTExd5vZs2fn9n3mNv/6c/jwYWOMMU899ZSRZPz8/Mw111xjBgwYYLp06WLq1KljHA6HefXVVwv0GAEAYIwxhHgAABSCs4V4mzdvzv0g+ssvv+RZ9+CDD+aGDH/9QJmVlWWGDx9uJJm4uDiTmZmZZ7u5c+ea/fv357uv5cuXm7CwMOPv72/27t2bZ92cOXNyA4Rly5blWffss8/mfng9lxBv9+7dpmTJkrkfart06WKef/55s3DhQnPixAmv23k8HtOyZUsjyQwfPtykpKTkrsvOzs4Nxtq2bZtnuy+//DI3vPh7CPLyyy973YczH+olmaZNm5ojR47krktKSjKlSpUykkzdunVNt27dTGpqau76X375xfj5+Rmn02l27dqV53bfeecdI8lUr17dbNiwIc+677//3pQsWdIEBASY33//Pc+6M71ERESYn3/+Oc+6J554wkgyNWvWzPe4nQnIli5dmm+dMca0bdvWSDLTpk3Lt+7EiRNmxYoVlttZSUlJMV9++WW+515WVpZ5+OGHjSTTpUsXrz1KMhMnTsyzbvHixbnh5p49e/Ksu++++4wkU7t2bbNv377c5ampqaZ79+65t/lvQrwDBw6Y0NBQ43A4zLp163LrzhbiDRw4MDfkPXjwYJ51r776qpFkatSoYXJycvKsmz59ujl9+nSeZR6Px4wfP95IMldccYXxeDx51v81RHrooYeM2+2+IP3ExcUZSblh718dPHgwX+h35jl4Lo+11X5YhXipqammcuXKRpIZMmSIMcaY5OTk3DD9v//9b57H5Zdffsn99/nOO+/kua1/CvEkmYYNG5oDBw7krktPTzfx8fFGkhkxYkS+/rw9D4z580uX4OBgU6JECbN169Z865OSksyWLVu8Pi4
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 37
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:37.808356Z",
|
||
|
|
"start_time": "2024-06-10T13:21:37.796358Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"cleaned_lines = [line for line in cleaned_lines if line not in closest_roads]\n",
|
||
|
|
"cleaned_lines.extend(new_segments)"
|
||
|
|
],
|
||
|
|
"id": "666f37cc24eb58e4",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": 38
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2024-06-10T13:21:39.663400Z",
|
||
|
|
"start_time": "2024-06-10T13:21:38.376155Z"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"cell_type": "code",
|
||
|
|
"source": [
|
||
|
|
"import heapq\n",
|
||
|
|
"G = nx.Graph()\n",
|
||
|
|
"\n",
|
||
|
|
"# Add edges to the graph from the cleaned lines\n",
|
||
|
|
"for line in cleaned_lines:\n",
|
||
|
|
" coords = list(line.coords)\n",
|
||
|
|
" for i in range(len(coords) - 1):\n",
|
||
|
|
" G.add_edge(coords[i], coords[i + 1], weight=Point(coords[i]).distance(Point(coords[i + 1])))\n",
|
||
|
|
"\n",
|
||
|
|
"# Function to find the paths between nearest points\n",
|
||
|
|
"def find_paths_between_nearest_points(G, nearest_points):\n",
|
||
|
|
" edges = []\n",
|
||
|
|
" for i, start_point in enumerate(nearest_points):\n",
|
||
|
|
" start = (start_point.x, start_point.y)\n",
|
||
|
|
" for end_point in nearest_points[i+1:]:\n",
|
||
|
|
" end = (end_point.x, end_point.y)\n",
|
||
|
|
" if nx.has_path(G, start, end):\n",
|
||
|
|
" path = nx.shortest_path(G, source=start, target=end, weight='weight')\n",
|
||
|
|
" path_edges = list(zip(path[:-1], path[1:]))\n",
|
||
|
|
" edges.extend((u, v, G[u][v]['weight']) for u, v in path_edges)\n",
|
||
|
|
" return edges\n",
|
||
|
|
"\n",
|
||
|
|
"# Find the edges used to connect the nearest points\n",
|
||
|
|
"edges = find_paths_between_nearest_points(G, nearest_points)\n",
|
||
|
|
"\n",
|
||
|
|
"# Create a graph from these edges\n",
|
||
|
|
"H = nx.Graph()\n",
|
||
|
|
"H.add_weighted_edges_from(edges)\n",
|
||
|
|
"\n",
|
||
|
|
"# Compute the Minimum Spanning Tree (MST) using Kruskal's algorithm\n",
|
||
|
|
"mst = nx.minimum_spanning_tree(H, weight='weight')\n",
|
||
|
|
"\n",
|
||
|
|
"# Perform pathfinding again on the MST to ensure shortest paths within the MST\n",
|
||
|
|
"final_edges = []\n",
|
||
|
|
"for u, v in mst.edges():\n",
|
||
|
|
" if nx.has_path(G, u, v):\n",
|
||
|
|
" path = nx.shortest_path(G, source=u, target=v, weight='weight')\n",
|
||
|
|
" path_edges = list(zip(path[:-1], path[1:]))\n",
|
||
|
|
" final_edges.extend((x, y, G[x][y]['weight']) for x, y in path_edges)\n",
|
||
|
|
"\n",
|
||
|
|
"# Create the final MST graph with these edges\n",
|
||
|
|
"final_mst = nx.Graph()\n",
|
||
|
|
"final_mst.add_weighted_edges_from(final_edges)\n",
|
||
|
|
"\n",
|
||
|
|
"# Identify and iteratively remove edges that do not have any nearest points and have one end with only one connection\n",
|
||
|
|
"nearest_points_tuples = [(point.x, point.y) for point in nearest_points]\n",
|
||
|
|
"\n",
|
||
|
|
"def find_edges_to_remove(graph):\n",
|
||
|
|
" edges_to_remove = []\n",
|
||
|
|
" for u, v in graph.edges():\n",
|
||
|
|
" if u not in nearest_points_tuples and v not in nearest_points_tuples:\n",
|
||
|
|
" if graph.degree(u) == 1 or graph.degree(v) == 1:\n",
|
||
|
|
" edges_to_remove.append((u, v))\n",
|
||
|
|
" return edges_to_remove\n",
|
||
|
|
"\n",
|
||
|
|
"edges_to_remove = find_edges_to_remove(final_mst)\n",
|
||
|
|
"\n",
|
||
|
|
"while edges_to_remove:\n",
|
||
|
|
" final_mst.remove_edges_from(edges_to_remove)\n",
|
||
|
|
" edges_to_remove = find_edges_to_remove(final_mst)\n",
|
||
|
|
"\n",
|
||
|
|
"# Function to generate a random color\n",
|
||
|
|
"def random_color():\n",
|
||
|
|
" return \"#\" + ''.join([random.choice('0123456789ABCDEF') for _ in range(6)])\n",
|
||
|
|
"\n",
|
||
|
|
"# Function to plot the graph with MST\n",
|
||
|
|
"# Function to plot the graph with MST\n",
|
||
|
|
"def plot_graph(G, mst, nearest_points):\n",
|
||
|
|
" fig, ax = plt.subplots(figsize=(40, 40)) # Increase the figure size\n",
|
||
|
|
"\n",
|
||
|
|
" pos = {node: node for node in G.nodes()}\n",
|
||
|
|
"\n",
|
||
|
|
" # Plot the graph edges\n",
|
||
|
|
" for (u, v) in G.edges():\n",
|
||
|
|
" x = [u[0], v[0]]\n",
|
||
|
|
" y = [u[1], v[1]]\n",
|
||
|
|
" ax.plot(x, y, color='blue', linewidth=1)\n",
|
||
|
|
"\n",
|
||
|
|
" # Plot the MST edges\n",
|
||
|
|
" for (u, v) in mst.edges():\n",
|
||
|
|
" x = [u[0], v[0]]\n",
|
||
|
|
" y = [u[1], v[1]]\n",
|
||
|
|
" ax.plot(x, y, color='green', linewidth=2) # Plot MST edges in green\n",
|
||
|
|
"\n",
|
||
|
|
" # Plot the nearest points\n",
|
||
|
|
" for point in nearest_points:\n",
|
||
|
|
" ax.plot(point.x, point.y, 'ro', markersize=8) # 'ro' for red circles, increase marker size\n",
|
||
|
|
"\n",
|
||
|
|
" ax.set_title('Graph with Incremental MST', fontsize=16)\n",
|
||
|
|
" ax.set_xlabel('Longitude', fontsize=14)\n",
|
||
|
|
" ax.set_ylabel('Latitude', fontsize=14)\n",
|
||
|
|
" plt.show()\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot the graph and MST\n",
|
||
|
|
"plot_graph(G, final_mst, nearest_points)"
|
||
|
|
],
|
||
|
|
"id": "4c7d8932a4fbb7c4",
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 4000x4000 with 1 Axes>"
|
||
|
|
],
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAADHkAAAxkCAYAAABvCGe/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeXicBb3/73fSUrakCFTWkIKgsrS0rCG0SGURKYWDgKJCKSIgCkhZiuz7jooHAXFBARFR4Ms5ekSQskkXJRIVkRQhCHKEhsWwHBogbfL74/m1IbZAgLTPtL3v68o1dJ5nZj4TUi+ZyWs+Vd3d3d0BAAAAAAAAAAAAAACgVNVlDwAAAAAAAAAAAAAAAIDIAwAAAAAAAAAAAAAAoCKIPAAAAAAAAAAAAAAAACqAyAMAAAAAAAAAAAAAAKACiDwAAAAAAAAAAAAAAAAqgMgDAAAAAAAAAAAAAACgAog8AAAAAAAAAAAAAAAAKoDIAwAAAAAAAAAAAAAAoAKIPAAAAAAAAAAAAAAAACqAyAMAAAAAgKXSPffck0MOOSQbb7xxVl555SyzzDJZddVVs/XWW+eII47I5MmT093dXfaYfTJmzJhUVVXlnnvuKXuU+TzxxBOpqqrKuuuu+65ve88996Sqqipjxozp15nWXXfdVFVV5eqrr+7X+6Xyzf13/8QTT/T5NmeccUaqqqpSVVWVD37wg+ns7HzLc5955pkMHDhw3vnXXXfdfOd0dXXl6quvzs4775zVVlstyyyzTFZZZZV85CMfyR577JGLLrpo3nxXX331vPt6N19+tgEAAAAAFl8Dyx4AAAAAAAAWpeeffz777bdffvOb3yRJ1l577YwaNSorrbRSXnrppTz00EO5/PLLc/nll2ezzTZLc3NzyRMvudZdd908+eST+fvf//6eIhCWTPfcc08+/vGPZ/vtt6+4cOn555/PL37xi+y9994LPH7NNddkzpw5b3n7V199NbvvvnvuvvvuJMnmm2+ej33sYxkwYEAef/zx3HbbbfnlL3+ZFVZYIUcccUQ22GCDTJgwYb77mTJlSlpbW7P++utn9OjR8x3fYIMN3uMzBAAAAACgbCIPAAAAAACWGi+++GJGjx6dRx55JBtuuGGuuOKKfPzjH5/vvIceeiiXXHJJbrjhhhKmXLKsvfbaaWlpyTLLLFP2KPC+bLnllvnDH/6QH/7wh28ZefzoRz/Ksssum49+9KN58MEH5zt+xhln5O67785aa62VX//619l00017HX/ppZdy8803Z80110ySjB49eoERx4EHHpjW1taMHj3a1g4AAAAAgCWMyAMAAAAAgKXGkUcemUceeSQf+tCHMm3atKy88soLPG/YsGG56qqr8qUvfWkRT7jkWWaZZbLhhhuWPQa8byNGjEhXV1duv/32PP3001lrrbV6Hb/vvvvyt7/9Lfvuu29mzpy5wPuYG46dfvrp8wUeSbLSSivloIMO6v/hAQAAAABYbFSXPQAAAAAAACwKra2tuf7665Mkl1xyyVsGHm+29dZbz3fdmDFjUlVVlXvuuSf33Xdfdt9993zwgx9MdXX1vE/Uf+WVV/L9738/e+21Vz784Q9nxRVXzIorrpjhw4fn5JNPzosvvrjAx1t33XVTVVWVJ554IrfccktGjx6dwYMHp7a2NmPGjMmtt976jjP/6U9/yl577ZUhQ4Zk2WWXzcYbb5xvfOMb6e7ufsfbznXppZemqqoqX/3qV+c7Nnbs2FRVVWWNNdaY7z6vvfbaVFVV5YADDph33RNPPJGqqqqsu+668667+uqrU1VVlSeffDJJst5666Wqqmre1z333DPf43Z2dubCCy/MJptskuWXXz6rrrpq9tprr7S0tPT5eb2TM844I1VVVTnjjDPy3HPP5fDDD88666yTQYMGZZ111smRRx75lv/ukuRvf/tbvvKVr+SjH/1oVlhhhQwePDgbb7xxvvKVr+Shhx6ad96bvydz5szJN7/5zWy22WapqalJVVXVfPf5pS99Keuvv36WW265rLTSSvnYxz6W6667boEzvPnn83e/+1122223rLrqqqmtrc3222+f++67b965t912W3bcccesvPLKqampyc4775zm5ua3fH7t7e05/fTTM3LkyNTW1maFFVbI8OHDc84552TWrFn98v0cM2bMvO069957b6+fizf/DD333HO59NJLM3bs2Ky33npZfvnlM3jw4Gy55Za58MIL89prr73l83g/DjrooMyZMyfXXHPNfMd++MMfzjvnrbS1tSVJVltttYUyHwAAAAAAiz+RBwAAAAAAS4X/+Z//SVdXV1ZeeeWMGzfufd/fjTfemDFjxuTxxx/PTjvtlJ133jnLLrtskuTPf/5zDj300EyZMiVrrLFGdt9994wePTrPPPNMzjvvvGy11VZ54YUX3vK+L7300uy11155/fXXM27cuGy88ca59957s9tuu+Xb3/72W97u9ttvT0NDQ2bMmJGdd945jY2N+dvf/pbjjjsuRx99dJ+f20477ZQkmTx5cq/rOzs789vf/jZJ8cvqf/nLX3odn3v+3Nu/lQ022CATJkzIiiuumCTZe++9M2HChHlfa6yxxnyPO3bs2Jx11lmpr6/PbrvtlhVXXDG33HJLtt122zzxxBN9fm598dRTT2XzzTfPzTffnK233jo777xzXnnllVx22WX5xCc+kc7Ozvluc/3112fTTTfNd77znbz22msZO3ZsdtpppwwaNChXXnllbrrppvlu093dnb322isnnnhiVl111eyxxx69tjvceOONGTFiRL73ve9l0KBBGTt2bLbccss0Nzdn/PjxbxsT/OpXv8p2222XZ555JjvvvHM22GCD/Pa3v83OO++cadOm5fLLL89uu+2W1157LZ/4xCey9tprZ/Lkydl+++3z2GOPzXd/Dz/8cEaMGJGzzjorzz77bEaPHp2ddtopzz33XE499dSMGjUqL7300vv+fn7yk5/MLrvskiRZffXVe/1c7LPPPvPOu/3223PUUUflwQcfzNChQ7Pnnntm6623ziOPPJITTjghO+ywQ15//fW3/P68V5///Oez3HLL5Uc/+lGv61955ZXceOONqa+vf9uf//r6+iTJlVdeuVDmAwAAAABg8Tew7AEAAAAAAGBReOCBB5Ikm2++eaqr3/9nIF1xxRW5/PLL85WvfGW+Y+uuu24mT56cj3/8470ea9asWfnyl7+ca6+9Nqeddlouv/zyBd73t771rVx33XXZb7/95l33s5/9LJ/73OdyzDHH5OMf/3iGDRs23+0uuOCCXHnllfnSl74077q77rorO+20Uy677LIcd9xxqaure8fntvHGG2ettdZKS0tLnn766ay11lpJkunTp+fVV1/NpptumgcffDCTJ0/uFSXceeedSd458hg9enRGjx6de+65J6+++mq+/vWv99rS8O+mTZuWzTbbLK2trfMCkNdeey177rlnbr/99px//vn57ne/+47Pq69++MMf5sADD8yVV145L9x56qmn0tjYmKamptx000353Oc+N+/8Bx54IAceeGBmz56dSy+9NIcffnivf+9PPvlknn/++fke5x//+Ee6urryl7/8JR/5yEd6HfvLX/6S8ePHp6qqKjfffHP22muvXve3++6750c/+lHGjBnTa3PKXN/4xjdy7bXXZv/995933bHHHptvfvObOeigg/LPf/4zv/nNb7LjjjsmSebMmZN99903N998cy688MJ8//vfn3e7jo6O7LHHHnnqqadyyimn5NRTT82gQYOSFD/TBx98cH7605/m6KOPnrfN4r1+P0844YRss802uf3227PhhhvO247z77bYYotMnz4922yzTa/r29vb89nPfja/+c1vcumll2bSpEkLvP17tfLKK+dTn/pUfvrTn+a+++7Ldtttl6T4+/nqq6/m2GOPfdv/fTniiCNy9NFH5/bbb8/QoUOzxx57ZJtttslmm22WTTfdNAMGDOjXeQEAAAAAWPzY5AEAAAAAwFJh7i/Zf/CDH1zg8T//+c858MAD5/uaMmXKAs/fYYcdFhh
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"execution_count": 39
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"metadata": {},
|
||
|
|
"cell_type": "code",
|
||
|
|
"outputs": [],
|
||
|
|
"execution_count": null,
|
||
|
|
"source": "",
|
||
|
|
"id": "e0b83c353f9d9191"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 2
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython2",
|
||
|
|
"version": "2.7.6"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 5
|
||
|
|
}
|