From addeddfe2a39338f59514b455b51e8ed8528d866 Mon Sep 17 00:00:00 2001 From: "Benjamin P. Stewart" Date: Wed, 30 Oct 2024 14:02:09 -0400 Subject: [PATCH] Updating friction surface generation code --- notebooks/TUTORIAL_Friction_Surface.ipynb | 6 +- notebooks/TUTORIAL_MCP_market_access.ipynb | 57 +++++++------ src/GOSTnetsraster/conversion_tables.py | 94 +++++++++++++++------- src/GOSTnetsraster/market_access.py | 42 +++++----- 4 files changed, 125 insertions(+), 74 deletions(-) diff --git a/notebooks/TUTORIAL_Friction_Surface.ipynb b/notebooks/TUTORIAL_Friction_Surface.ipynb index 3f8d712..f015027 100644 --- a/notebooks/TUTORIAL_Friction_Surface.ipynb +++ b/notebooks/TUTORIAL_Friction_Surface.ipynb @@ -172,9 +172,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python (geog)", + "display_name": "gn", "language": "python", - "name": "geog" + "name": "gn" }, "language_info": { "codemirror_mode": { @@ -186,7 +186,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/notebooks/TUTORIAL_MCP_market_access.ipynb b/notebooks/TUTORIAL_MCP_market_access.ipynb index dfa813d..00107da 100644 --- a/notebooks/TUTORIAL_MCP_market_access.ipynb +++ b/notebooks/TUTORIAL_MCP_market_access.ipynb @@ -15,7 +15,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -35,38 +35,42 @@ "from scipy.ndimage import generic_filter\n", "from pandana.loaders import osm\n", "\n", + "sys.path.insert(0, r\"C:\\WBG\\Work\\Code\\GOSTrocks\\src\")\n", "sys.path.append(\"../src\")\n", "\n", - "import GOSTNetsRaster.market_access as ma" + "import GOSTnetsraster.market_access as ma" ] }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { - "ename": "DriverError", - "evalue": "../tutorial_data/destinations.shp: No such file or directory", + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\pyogrio\\core.py:35: RuntimeWarning: Could not detect GDAL data files. Set GDAL_DATA environment variable to the correct path.\n", + " _init_gdal_data()\n" + ] + }, + { + "ename": "DataSourceError", + "evalue": "../tutorial_data\\destinations.shp: No such file or directory", "output_type": "error", "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mCPLE_OpenFailedError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32mfiona/_shim.pyx\u001b[0m in \u001b[0;36mfiona._shim.gdal_open_vector\u001b[0;34m()\u001b[0m\n", - "\u001b[0;32mfiona/_err.pyx\u001b[0m in \u001b[0;36mfiona._err.exc_wrap_pointer\u001b[0;34m()\u001b[0m\n", - "\u001b[0;31mCPLE_OpenFailedError\u001b[0m: ../tutorial_data/destinations.shp: No such file or directory", - "\nDuring handling of the above exception, another exception occurred:\n", - "\u001b[0;31mDriverError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m/tmp/ipykernel_10421/1075862729.py\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0minD\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdests\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7\u001b[0m \u001b[0minR\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrasterio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfriction_surface\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mfrictionD\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minR\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/geog/lib/python3.9/site-packages/geopandas/io/file.py\u001b[0m in \u001b[0;36m_read_file\u001b[0;34m(filename, bbox, mask, rows, **kwargs)\u001b[0m\n\u001b[1;32m 199\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 200\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mfiona_env\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 201\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mreader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_or_bytes\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mfeatures\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 202\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 203\u001b[0m \u001b[0;31m# In a future Fiona release the crs attribute of features will\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/geog/lib/python3.9/site-packages/fiona/env.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m 406\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 407\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlocal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_env\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 408\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 409\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 410\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/geog/lib/python3.9/site-packages/fiona/__init__.py\u001b[0m in \u001b[0;36mopen\u001b[0;34m(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, **kwargs)\u001b[0m\n\u001b[1;32m 254\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 255\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mmode\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 256\u001b[0;31m c = Collection(path, mode, driver=driver, encoding=encoding,\n\u001b[0m\u001b[1;32m 257\u001b[0m layer=layer, enabled_drivers=enabled_drivers, **kwargs)\n\u001b[1;32m 258\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mmode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'w'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/.conda/envs/geog/lib/python3.9/site-packages/fiona/collection.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, **kwargs)\u001b[0m\n\u001b[1;32m 160\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 161\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msession\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mSession\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 162\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msession\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 163\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmode\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m'a'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'w'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 164\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msession\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mWritingSession\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32mfiona/ogrext.pyx\u001b[0m in \u001b[0;36mfiona.ogrext.Session.start\u001b[0;34m()\u001b[0m\n", - "\u001b[0;32mfiona/_shim.pyx\u001b[0m in \u001b[0;36mfiona._shim.gdal_open_vector\u001b[0;34m()\u001b[0m\n", - "\u001b[0;31mDriverError\u001b[0m: ../tutorial_data/destinations.shp: No such file or directory" + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mDataSourceError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn[4], line 6\u001b[0m\n\u001b[0;32m 2\u001b[0m dests \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(tutorial_folder, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdestinations.shp\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 3\u001b[0m friction_surface \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mpath\u001b[38;5;241m.\u001b[39mjoin(tutorial_folder, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mglobal_friction_surface.tif\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m----> 6\u001b[0m inD \u001b[38;5;241m=\u001b[39m gpd\u001b[38;5;241m.\u001b[39mread_file(dests)\n\u001b[0;32m 7\u001b[0m inR \u001b[38;5;241m=\u001b[39m rasterio\u001b[38;5;241m.\u001b[39mopen(friction_surface)\n\u001b[0;32m 8\u001b[0m frictionD \u001b[38;5;241m=\u001b[39m inR\u001b[38;5;241m.\u001b[39mread()[\u001b[38;5;241m0\u001b[39m,:,:]\n", + "File \u001b[1;32mC:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\geopandas\\io\\file.py:294\u001b[0m, in \u001b[0;36m_read_file\u001b[1;34m(filename, bbox, mask, columns, rows, engine, **kwargs)\u001b[0m\n\u001b[0;32m 291\u001b[0m from_bytes \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[0;32m 293\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpyogrio\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m--> 294\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _read_file_pyogrio(\n\u001b[0;32m 295\u001b[0m filename, bbox\u001b[38;5;241m=\u001b[39mbbox, mask\u001b[38;5;241m=\u001b[39mmask, columns\u001b[38;5;241m=\u001b[39mcolumns, rows\u001b[38;5;241m=\u001b[39mrows, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs\n\u001b[0;32m 296\u001b[0m )\n\u001b[0;32m 298\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m engine \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfiona\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m 299\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m pd\u001b[38;5;241m.\u001b[39mapi\u001b[38;5;241m.\u001b[39mtypes\u001b[38;5;241m.\u001b[39mis_file_like(filename):\n", + "File \u001b[1;32mC:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\geopandas\\io\\file.py:547\u001b[0m, in \u001b[0;36m_read_file_pyogrio\u001b[1;34m(path_or_bytes, bbox, mask, rows, **kwargs)\u001b[0m\n\u001b[0;32m 538\u001b[0m warnings\u001b[38;5;241m.\u001b[39mwarn(\n\u001b[0;32m 539\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124minclude_fields\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m and \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mignore_fields\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m keywords are deprecated, and \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 540\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mwill be removed in a future release. You can use the \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mcolumns\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m keyword \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 543\u001b[0m stacklevel\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m3\u001b[39m,\n\u001b[0;32m 544\u001b[0m )\n\u001b[0;32m 545\u001b[0m kwargs[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcolumns\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m kwargs\u001b[38;5;241m.\u001b[39mpop(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124minclude_fields\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m--> 547\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m pyogrio\u001b[38;5;241m.\u001b[39mread_dataframe(path_or_bytes, bbox\u001b[38;5;241m=\u001b[39mbbox, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n", + "File \u001b[1;32mC:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\pyogrio\\geopandas.py:265\u001b[0m, in \u001b[0;36mread_dataframe\u001b[1;34m(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, fid_as_index, use_arrow, on_invalid, arrow_to_pandas_kwargs, **kwargs)\u001b[0m\n\u001b[0;32m 260\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m use_arrow:\n\u001b[0;32m 261\u001b[0m \u001b[38;5;66;03m# For arrow, datetimes are read as is.\u001b[39;00m\n\u001b[0;32m 262\u001b[0m \u001b[38;5;66;03m# For numpy IO, datetimes are read as string values to preserve timezone info\u001b[39;00m\n\u001b[0;32m 263\u001b[0m \u001b[38;5;66;03m# as numpy does not directly support timezones.\u001b[39;00m\n\u001b[0;32m 264\u001b[0m kwargs[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdatetime_as_string\u001b[39m\u001b[38;5;124m\"\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m--> 265\u001b[0m result \u001b[38;5;241m=\u001b[39m read_func(\n\u001b[0;32m 266\u001b[0m path_or_buffer,\n\u001b[0;32m 267\u001b[0m layer\u001b[38;5;241m=\u001b[39mlayer,\n\u001b[0;32m 268\u001b[0m encoding\u001b[38;5;241m=\u001b[39mencoding,\n\u001b[0;32m 269\u001b[0m columns\u001b[38;5;241m=\u001b[39mcolumns,\n\u001b[0;32m 270\u001b[0m read_geometry\u001b[38;5;241m=\u001b[39mread_geometry,\n\u001b[0;32m 271\u001b[0m force_2d\u001b[38;5;241m=\u001b[39mgdal_force_2d,\n\u001b[0;32m 272\u001b[0m skip_features\u001b[38;5;241m=\u001b[39mskip_features,\n\u001b[0;32m 273\u001b[0m max_features\u001b[38;5;241m=\u001b[39mmax_features,\n\u001b[0;32m 274\u001b[0m where\u001b[38;5;241m=\u001b[39mwhere,\n\u001b[0;32m 275\u001b[0m bbox\u001b[38;5;241m=\u001b[39mbbox,\n\u001b[0;32m 276\u001b[0m mask\u001b[38;5;241m=\u001b[39mmask,\n\u001b[0;32m 277\u001b[0m fids\u001b[38;5;241m=\u001b[39mfids,\n\u001b[0;32m 278\u001b[0m sql\u001b[38;5;241m=\u001b[39msql,\n\u001b[0;32m 279\u001b[0m sql_dialect\u001b[38;5;241m=\u001b[39msql_dialect,\n\u001b[0;32m 280\u001b[0m return_fids\u001b[38;5;241m=\u001b[39mfid_as_index,\n\u001b[0;32m 281\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs,\n\u001b[0;32m 282\u001b[0m )\n\u001b[0;32m 284\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m use_arrow:\n\u001b[0;32m 285\u001b[0m meta, table \u001b[38;5;241m=\u001b[39m result\n", + "File \u001b[1;32mC:\\wbg\\Anaconda3\\envs\\gn\\Lib\\site-packages\\pyogrio\\raw.py:198\u001b[0m, in \u001b[0;36mread\u001b[1;34m(path_or_buffer, layer, encoding, columns, read_geometry, force_2d, skip_features, max_features, where, bbox, mask, fids, sql, sql_dialect, return_fids, datetime_as_string, **kwargs)\u001b[0m\n\u001b[0;32m 59\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Read OGR data source into numpy arrays.\u001b[39;00m\n\u001b[0;32m 60\u001b[0m \n\u001b[0;32m 61\u001b[0m \u001b[38;5;124;03mIMPORTANT: non-linear geometry types (e.g., MultiSurface) are converted\u001b[39;00m\n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 194\u001b[0m \n\u001b[0;32m 195\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 196\u001b[0m dataset_kwargs \u001b[38;5;241m=\u001b[39m _preprocess_options_key_value(kwargs) \u001b[38;5;28;01mif\u001b[39;00m kwargs \u001b[38;5;28;01melse\u001b[39;00m {}\n\u001b[1;32m--> 198\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m ogr_read(\n\u001b[0;32m 199\u001b[0m get_vsi_path_or_buffer(path_or_buffer),\n\u001b[0;32m 200\u001b[0m layer\u001b[38;5;241m=\u001b[39mlayer,\n\u001b[0;32m 201\u001b[0m encoding\u001b[38;5;241m=\u001b[39mencoding,\n\u001b[0;32m 202\u001b[0m columns\u001b[38;5;241m=\u001b[39mcolumns,\n\u001b[0;32m 203\u001b[0m read_geometry\u001b[38;5;241m=\u001b[39mread_geometry,\n\u001b[0;32m 204\u001b[0m force_2d\u001b[38;5;241m=\u001b[39mforce_2d,\n\u001b[0;32m 205\u001b[0m skip_features\u001b[38;5;241m=\u001b[39mskip_features,\n\u001b[0;32m 206\u001b[0m max_features\u001b[38;5;241m=\u001b[39mmax_features \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;241m0\u001b[39m,\n\u001b[0;32m 207\u001b[0m where\u001b[38;5;241m=\u001b[39mwhere,\n\u001b[0;32m 208\u001b[0m bbox\u001b[38;5;241m=\u001b[39mbbox,\n\u001b[0;32m 209\u001b[0m mask\u001b[38;5;241m=\u001b[39m_mask_to_wkb(mask),\n\u001b[0;32m 210\u001b[0m fids\u001b[38;5;241m=\u001b[39mfids,\n\u001b[0;32m 211\u001b[0m sql\u001b[38;5;241m=\u001b[39msql,\n\u001b[0;32m 212\u001b[0m sql_dialect\u001b[38;5;241m=\u001b[39msql_dialect,\n\u001b[0;32m 213\u001b[0m return_fids\u001b[38;5;241m=\u001b[39mreturn_fids,\n\u001b[0;32m 214\u001b[0m dataset_kwargs\u001b[38;5;241m=\u001b[39mdataset_kwargs,\n\u001b[0;32m 215\u001b[0m datetime_as_string\u001b[38;5;241m=\u001b[39mdatetime_as_string,\n\u001b[0;32m 216\u001b[0m )\n", + "File \u001b[1;32mpyogrio\\\\_io.pyx:1240\u001b[0m, in \u001b[0;36mpyogrio._io.ogr_read\u001b[1;34m()\u001b[0m\n", + "File \u001b[1;32mpyogrio\\\\_io.pyx:220\u001b[0m, in \u001b[0;36mpyogrio._io.ogr_open\u001b[1;34m()\u001b[0m\n", + "\u001b[1;31mDataSourceError\u001b[0m: ../tutorial_data\\destinations.shp: No such file or directory" ] } ], @@ -147,6 +151,13 @@ "ma.generate_market_sheds(inR, inD, out_file=outfile)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, @@ -155,9 +166,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python (geog)", + "display_name": "gn", "language": "python", - "name": "geog" + "name": "gn" }, "language_info": { "codemirror_mode": { @@ -169,7 +180,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/src/GOSTnetsraster/conversion_tables.py b/src/GOSTnetsraster/conversion_tables.py index 145d864..98d712d 100755 --- a/src/GOSTnetsraster/conversion_tables.py +++ b/src/GOSTnetsraster/conversion_tables.py @@ -14,33 +14,46 @@ 'living_street':10, 'service':10 } - -copernicus_landcover = { - 0:0.1, - 111: 3.24, - 113: 3.24, - 112: 1.62, - 114: 4.00, - 115: 3.24, - 116: 3.24, - 121: 3.75, - 123: 3.75, - 122: 2.10, - 124: 4.00, - 125: 3.75, - 126: 3.75, - 20: 3.00, - 30: 4.20, - 90: 2.00, - 100: 4.86, - 60: 4.86, - 40: 2.50, - 50: 5, - 70: 1.62, - 80: 0.5, - 200: 0, +esaacci_landcover = { + 0: 0.1, #No data + 10: 2.50, #Cropland, rainfed + 11: 4.00, # ???Herbaceous cover + 12: 4.20, # ???Tree or shrub cover + 20: 3.24, #Cropland, irrigated or post-flooding + 30: 2.00, #Mosaic cropland (>50%) / natural vegetation (tree, shrub, herbaceous cover) (<50%) + 40: 3.24, #Mosaic natural vegetation (tree, shrub, herbaceous cover) (>50%) / cropland (<50%) + 50: 4.00, #Tree cover, broadleaved, evergreen, closed to open (>15%) + 60: 3.50, #Tree cover, broadleaved, deciduous, closed to open (>15%) + 61: 3.00, #Tree cover, broadleaved, deciduous, closed (>40%) + 62: 2.50, #Tree cover, broadleaved, deciduous, open (15-40%) + 70: 3.50, #Tree cover, needleleaved, evergreen, closed to open (>15%) + 71: 3.00, #Tree cover, needleleaved, evergreen, closed (>40%) + 72: 3.24, #Tree cover, needleleaved, evergreen, open (15-40%) + 80: 3.50, #Tree cover, needleleaved, deciduous, closed to open (>15%) + 81: 3.00, #Tree cover, needleleaved, deciduous, closed (>40%) + 82: 3.24, #Tree cover, needleleaved, deciduous, open (15-40%) + 90: 3.24, #Tree cover, mixed leaf type (broadleaved and needleleaved) + 100: 3.00, #Mosaic tree and shrub (>50%) / herbaceous cover (<50%) + 110: 3.00, #Mosaic herbaceous cover (>50%) / tree and shrub (<50%) + 120: 4.20, #Shrubland + 121: 3.24, #Shrubland evergreen + 122: 3.24, #Shrubland deciduous + 130: 4.86, #Grassland + 140: 0.00, # ???Lichens and mosses + 150: 4.20, #Sparse vegetation (tree, shrub, herbaceous cover) (<15%) + 151: 4.00, #Sparse tree (<15%) + 152: 4.00, #Sparse shrub (<15%) + 153: 4.00, #Sparse herbaceous cover (<15%) + 160: 2.00, #Tree cover, flooded, fresh or brakish water + 170: 2.00, #Tree cover, flooded, saline water + 180: 2.00, #Shrub or herbaceous cover, flooded, fresh/saline/brakish water + 190: 5.00, #Urban areas + 200: 4.75, #Bare areas + 201: 4.50, #Consolidated bare areas + 202: 4.50, #Unconsolidated bare areas + 210: 0.50, #Water bodies + 220: 1.00, #Permanent snow and ice } - ''' https://www.nature.com/articles/s41591-020-1059-1#MOESM2 evergreen needleleaf forest = 3.24 evergreen broadleaf forest = 1.62 @@ -57,4 +70,31 @@ cropland/natural vegetation = 3.24 snow and ice = 1.62 barren or sparsely vegetated = 3.00 -''' \ No newline at end of file +''' +copernicus_landcover = { + 0:0.1, + 10: 2.50, + 20: 3.00, + 30: 4.20, + 40: 2.50, + 50: 5, + 60: 4.86, + 70: 1.62, + 80: 0.5, + 90: 2.00, + 100: 4.86, + 111: 3.24, + 113: 3.24, + 112: 1.62, + 114: 4.00, + 115: 3.24, + 116: 3.24, + 121: 3.75, + 123: 3.75, + 122: 2.10, + 124: 4.00, + 125: 3.75, + 126: 3.75, + 200: 0, +} + diff --git a/src/GOSTnetsraster/market_access.py b/src/GOSTnetsraster/market_access.py index fbae74d..24df6dc 100644 --- a/src/GOSTnetsraster/market_access.py +++ b/src/GOSTnetsraster/market_access.py @@ -12,7 +12,7 @@ import osmnx as ox import GOSTnets as gn import skimage.graph as graph -import GOSTRocks.rasterMisc as rMisc +import GOSTrocks.rasterMisc as rMisc from skimage.graph import _mcp from rasterio.mask import mask @@ -105,7 +105,9 @@ def name_mcp_dests(inH, destinations): destinations.loc[idx, 'MCP_DESTS_NAME'] = "_".join([str(xx) for xx in inH.index(x.x, x.y)]) return(destinations) -def generate_roads_lc_friction(lc_file, sel_roads, lc_travel_table=None, min_lc_val=0.01, min_road_speed=0.01, speed_col='speed', resolution=100, out_file = ''): +def generate_roads_lc_friction(lc_file, sel_roads, lc_travel_table=None, min_lc_val=0.01, + min_road_speed=0.01, speed_col='speed', + resolution=100, out_file = ''): ''' Combine a landcover dataset and a road network dataset to create a friction surface. See generate_lc_friction and generate_road_friction for more details @@ -116,17 +118,17 @@ def generate_roads_lc_friction(lc_file, sel_roads, lc_travel_table=None, min_lc_ if not lc_travel_table: lc_travel_table = speed_tables.copernicus_landcover lc_friction = generate_lc_friction(lc_file, lc_travel_table = lc_travel_table, min_val = min_lc_val, resolution = resolution) - road_friction = generate_network_raster(lc_file, sel_roads, min_speed=min_road_speed, speed_col=speed_col, resolution=resolution) + road_friction = generate_road_friction(lc_file, sel_roads, no_road_speed=min_road_speed, speed_col=speed_col, resolution=resolution) #Stack frictions and find minimum stacked_friction = np.dstack([road_friction, lc_friction[0,:,:]]) combo_friction = np.amin(stacked_friction, axis=2) - combo_friction[combo_friction == inf] = 65000 + combo_friction[combo_friction == inf] = None + combo_friction = combo_friction.astype('float') out_meta = lc_file.meta.copy() - combo_friction = combo_friction.astype('uint16') - out_meta['dtype'] = combo_friction.dtype - + out_meta.update(dtype=combo_friction.dtype) + if out_file != '': with rasterio.open(out_file, 'w', **out_meta) as outR: outR.write_band(1, combo_friction) @@ -150,32 +152,30 @@ def generate_lc_friction(lc_file, lc_travel_table=None, min_val=0.01, resolution lc_data = lc_file.read() res = np.vectorize(lc_travel_table.get)(lc_data) res = res.astype(float) - res = np.nan_to_num(res, nan=0.01) - res = resolution / (res * 1000 / (60 * 60)) # km/h --> m/s * resolution of image in metres + res = np.nan_to_num(res, nan=min_val) + res = 1/((res * 1000)/60) # km/h --> m/min --> minutes/m return(res) -def generate_road_friction(inH, sel_roads, min_speed=0.01, speed_col='speed', resolution=100): +def generate_road_friction(inH, sel_roads, no_road_speed=0.01, speed_col='speed', resolution=100): ''' Create raster with network travel times from a road network that measures seconds to cross a cell INPUTS inH [rasterio object] - template raster used to define raster shape, resolution, crs, etc. sel_roads [geopandas dataframe] - road network to burn into raster - [optional] min_speed [int] - minimum travel speed for areas without roads + [optional] no_road_speed [int] - minimum travel speed for areas without roads [optional] speed_col [string] - column in sel_roads that defines the speed in KM/h [optional] resolution [int] - resolution of the raster in metres RETURNS [numpy array] ''' - # create a copy of inH with value set to slowest walking speed - distance_data = np.zeros(inH.shape) # burn the speeds into the distance_data using the road network sel_roads = sel_roads.sort_values([speed_col]) shapes = ((row['geometry'], row[speed_col]) for idx, row in sel_roads.iterrows()) - speed_image = features.rasterize(shapes, out_shape=inH.shape, transform=inH.transform, fill=min_speed) + res = features.rasterize(shapes, out_shape=inH.shape, transform=inH.transform, fill=no_road_speed) # convert to a version that claculates the seconds to cross each cell - traversal_time = resolution / (speed_image * 1000 / (60 * 60)) # km/h --> m/s * resolution of image in metres - return(traversal_time) + res = 1/((res * 1000)/60) # km/h --> m/min --> minutes/m + return(res) def calculate_travel_time(inH, mcp, destinations, out_raster = ''): ''' Calculate travel time raster @@ -272,10 +272,10 @@ def generate_feature_vectors(network_r, mcp, inH, threshold, featIdx='tempID', v gs = gpd.GeoSeries()#Create an empty GeoSeries. for x in all_shapes: - gsTemp = gpd.GeoSeries(x) - gs = gs.append(gsTemp) - - gs = gs.reset_index(drop=True)#Reset the index of 'gs' + gsTemp = gpd.GeoSeries(x) + gs = gs.append(gsTemp) + + gs = gs.reset_index(drop=True) # Reset the index of 'gs' # Geometry validity check and correction @@ -290,7 +290,7 @@ def generate_feature_vectors(network_r, mcp, inH, threshold, featIdx='tempID', v # This code block can be omitted if it's redundant. for i, geom in enumerate(gs): print('Geometry No. {}: Validity = {} | THS = {} minutes | Poly Count = {}'.format(i, geom.is_valid, thresh, polyCount)) - + union = gs.unary_union complete_shapes.append([union, thresh, row[featIdx]])