diff --git a/README.md b/README.md index 013282f1..2b71db91 100644 --- a/README.md +++ b/README.md @@ -131,6 +131,8 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e #### Supported Spatial Functions with Supported Inputs +βœ… = native backend    πŸ”„ = accepted (CPU fallback) + ------- ### **Classification** @@ -183,7 +185,7 @@ In the GIS world, rasters are used for representing continuous phenomena (e.g. e | Name | Description | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray | |:----------:|:------------|:----------------------:|:--------------------:|:-------------------:|:------:| -| [A* Pathfinding](xrspatial/pathfinding.py) | Finds the least-cost path between two cells on a cost surface | βœ…οΈ | | | | +| [A* Pathfinding](xrspatial/pathfinding.py) | Finds the least-cost path between two cells on a cost surface | βœ…οΈ | βœ… | πŸ”„ | πŸ”„ | ---------- diff --git a/docs/source/user_guide/pathfinding.ipynb b/docs/source/user_guide/pathfinding.ipynb index 1bdb7da5..0e42ae77 100644 --- a/docs/source/user_guide/pathfinding.ipynb +++ b/docs/source/user_guide/pathfinding.ipynb @@ -3,16 +3,7 @@ { "cell_type": "markdown", "metadata": {}, - "source": [ - "## Pathfinding" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Xarray-spatial's Pathfinding provides a comprehensive tool for finding the shortest path from one point to another in a raster that can contain any level of complex boundaries or obstacles amidst an interconnected set of traversable path segments." - ] + "source": "## Pathfinding\n\nXarray-spatial provides A\\* pathfinding for finding optimal routes through raster surfaces.\nPaths can be computed using geometric distance alone (unweighted) or weighted by a\nfriction/cost surface, which makes the algorithm find the *least-cost* path rather\nthan the geometrically shortest one.\n\nAll four backends are supported:\n\n| Backend | Strategy |\n|---------|----------|\n| **NumPy** | Numba-jitted kernel (fast, in-memory) |\n| **Dask** | Sparse Python A\\* with LRU chunk cache β€” loads chunks on demand so the full grid never needs to fit in RAM |\n| **CuPy** | CPU fallback (transfers to numpy, runs the numba kernel, transfers back) |\n| **Dask+CuPy** | Same sparse A\\* as Dask, with automatic cupy-to-numpy chunk conversion |\n\n**Note:** `snap_start` and `snap_goal` are not supported with Dask-backed arrays\nbecause the brute-force nearest-pixel scan is O(h*w). Ensure the start and goal\npixels are valid before calling `a_star_search`.\n\n- [Unweighted A\\*](#Unweighted-A*): find the shortest geometric path through a line network\n- [Weighted A\\*](#Weighted-A*): find the least-cost path through a friction surface\n- [Dask (out-of-core)](#Dask-(out-of-core)): pathfinding on datasets that don't fit in RAM" }, { "cell_type": "markdown", @@ -23,118 +14,44 @@ }, { "cell_type": "code", - "execution_count": 1, - "metadata": { - "execution": { - "iopub.execute_input": "2026-01-23T01:07:47.029270Z", - "iopub.status.busy": "2026-01-23T01:07:47.029151Z", - "iopub.status.idle": "2026-01-23T01:07:47.924441Z", - "shell.execute_reply": "2026-01-23T01:07:47.923802Z" - } - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "import pandas as pd\n", - "\n", - "import datashader as ds\n", - "from datashader.transfer_functions import shade\n", - "from datashader.transfer_functions import stack\n", - "from datashader.transfer_functions import dynspread\n", - "from datashader.transfer_functions import set_background\n", - "from datashader.colors import Elevation\n", - "\n", - "import xrspatial" - ] - }, - { - "cell_type": "markdown", + "execution_count": null, "metadata": {}, - "source": [ - "To download the examples data, run the command `xrspatial examples` in your terminal. All the data will be stored in your current directory inside a folder named `xrspatial-examples`." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### A*" - ] + "outputs": [], + "source": "import numpy as np\nimport pandas as pd\nimport xarray as xr\n\nimport datashader as ds\nfrom datashader.transfer_functions import shade\nfrom datashader.transfer_functions import stack\nfrom datashader.transfer_functions import dynspread\nfrom datashader.transfer_functions import set_background\nfrom datashader.colors import Elevation\n\nimport xrspatial\nfrom xrspatial import a_star_search, generate_terrain, slope, cost_distance" }, { "cell_type": "markdown", "metadata": {}, "source": [ - "A* is an informed search algorithm, or a best-first search, meaning that it is formulated in terms of weighted graphs: starting from a specific starting node of a graph, it aims to find a path to the given goal node having the smallest cost (min distance travelled, shortest time, ...). \n", + "### Unweighted A*\n", + "\n", + "A\\* is an informed search algorithm that finds the shortest path from a start\n", + "to a goal through a graph. In the unweighted case, edge costs are purely\n", + "geometric (Euclidean distance scaled by cellsize).\n", "\n", - "The `xrspatial.a_star_search` function calculates the shortest path in pixel space from a start location to a goal location through a given aggregate surface graph. The graph should be a line raster which contains crossable and non-crossable (a.k.a walls or barrieres) values. Note that both start and goal are in (lon, lat), or (x, y) coordinate space and must be within the graph. `xrspatial.a_star_search` provides 2 separate options, `snap_start` and `snap_goal`, which can be set to true to snap locations to the nearest valid value before beginning pathfinding. It also provides a `connectivity` option to indicate neighborhood structure. This value can be set to 4 or 8 for 4-connectivity or 8-connectivity." + "We'll demonstrate this on a synthetic line raster with barriers." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Let's generate a fake line raster and find the shortest path with A*.\n", - "\n", - "- First, we'll generate a line raster by setting up a pandas DataFrame specifying the line coordinates.\n", - "- Then, we'll aggregate that into a lines raster with Canvas.line\n", - "- Once we have that, we'll choose a start and goal point to put into the a* pathfinding function.\n", - "- For visualization, we'll also aggregate those points and render them in an image together with the lines." + "#### Create a line raster" ] }, { "cell_type": "code", - "execution_count": 2, - "metadata": { - "execution": { - "iopub.execute_input": "2026-01-23T01:07:47.926119Z", - "iopub.status.busy": "2026-01-23T01:07:47.925813Z", - "iopub.status.idle": "2026-01-23T01:07:49.977084Z", - "shell.execute_reply": "2026-01-23T01:07:49.976465Z" - } - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAA3CElEQVR4Ae3da5Iju3EG0JHDoQV7A9YKvGD7h6ykLntyqslukI1CAchTEY4Lsot4nIRC/gTI/tuvX7/++a//8RAgQIAAAQIECBAgQOB0gf84fQQDECBAgAABAgQIECBA4C8BAcRWIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECAggNgDBAgQIECAAAECBAgMExBAhlEbiAABAgQIECBAgAABAcQeIECAAAECBAgQIEBgmIAAMozaQAQIECBAgAABAgQICCD2AAECBAgQIECAAAECwwQEkGHUBiJAgAABAgQIECBAQACxBwgQIECAAAECBAgQGCYggAyjNhABAgQIECBAgAABAgKIPUCAAAECBAgQIECAwDABAWQYtYEIECBAgAABAgQIEBBA7AECBAgQIECAAAECBIYJCCDDqA1EgAABAgQIECBAgIAAYg8QIECAAAECBAgQIDBMQAAZRm0gAgQIECBAgAABAgQEEHuAAAECBAgQIECAAIFhAgLIMGoDESBAgAABAgQIECDwnwgIEPgs8M/PX/3624PvfEWAAAECBAgQIPCagBOQ17y8XUDgUfiIZT/7vgCJJRIgQIAAAQIEugkIIN0odbSDwHch47u/72BgDQQIECBAgACBMwUEkDN19b2UQGu4aH1vqcWbLAECBAgQIEBgkIAAMgjaMAQIECBAgAABAgQI/Prlv4RuFxB4Q+B///Ffb/zKTwgQIECgksDf//E/lZZrrQSaBQSQZiovEvgt4N9Uflto1RWIIO5fC3Xrb+WPBfwHVI9dfEsgC7iClTW0CRAgQIAAAQJvCtxDuWD+JqCflREQQMqU2kK/E2j9//PR+t534/k7AQIECOwjcA8f+6zISgicJyCAnGer5wUFvgsX8XfH6wsW1pQJECBwkkD8e4LwcRKubrcV8N8B2ba0FvauQISMR/+ndu/hJI7W7yHEMfu7yn5HgACB9QUEj/VraAXXCDgBucbdqJMLRNiI//m/f/0nW/d2nnIEjxxE8t+0CRAgQGB/AeFj/xpb4XkCAsh5tnouICCEFCiyJRIgQCAJRPAQPhKIJoE3BFzBegPNTwhkgRxCou0hQIAAgT0FBI8962pV4wWcgIw3N+KGAhE8chDZcImWRIAAgdICwkfp8lt8ZwEBpDOo7moLCCG162/1BAjsJ+DK1X41taLrBVzBur4GZrCZQA4h0fYQIECAwJoCTj3WrJtZzy/gBGT+GpnhggIRPHIQWXAJpkyAAIHSAsJH6fJb/MkCAsjJwLqvLSCE1K6/1RMgsJ6AK1fr1cyM1xNwBWu9mpnxYgI5hETbQ4AAAQJzCjj1mLMuZrWfgBOQ/WpqRRMKRPDIQWTCKZoSAQIESgsIH6XLb/GDBQSQweCGqy0ghNSuv9UTIDCfgCtX89XEjPYXcAVr/xpb4WQCOYRE20OAAAEC1wg49bjG3agEnIDYAwQuEIjgkYPIBVMwJAECBEoLCB+ly2/xFwsIIBcXwPC1BYSQ2vW3egIExgu4cjXe3IgEjgKuYB1FfCYwWCCHkGh7CBAgQOAcAace57jqlcCrAk5AXhXzPoETBCJ45CBywhC6JECAQGkB4aN0+S1+MgEBZLKCmE5tASGkdv2tngCB/gKuXPU31SOBnwq4gvVTQb8n0Fkgh5BoewgQIEDgPQGnHu+5+RWBswWcgJwtrH8CbwhE8MhB5I0u/IQAAQKlBYSP0uW3+MkFBJDJC2R6tQWEkNr1t3oCBF4XcOXqdTO/IDBawBWs0eLGI/CiQA4h0fYQIECAwGMBpx6PXXxLYDYBJyCzVcR8CDwQiOCRg8iDV3xFgACB0gLCR+nyW/xiAgLIYgUz3doCQkjt+ls9AQKfBVy5+mziGwKzC7iCNXuFzI/AQSCHkGh7CBAgUFXAqUfVylv36gJOQFavoPmXFIjgkYNISQSLJkCgtIDwUbr8Fr+4gACyeAFNv7aAEFK7/lZPoKKAK1cVq27Nuwm4grVbRa2nnEAOIdH2ECBAYFcBpx67Vta6qgk4AalWcevdUiCCRw4iWy7SoggQKC0gfJQuv8VvJiCAbFZQy6ktIITUrr/VE9hRwJWrHatqTdUFXMGqvgOsfzuBHEKi7SFAgMCqAk49Vq2ceRP4WsAJyNc+/kpgSYEIHjmILLkIkyZAoLSA8FG6/Ba/uYAAsnmBLa+2gBBSu/5WT2BFAVeuVqyaORN4TcAVrNe8vE1gOYEcQqLtIUCAwKwCTj1mrYx5Eegr4ASkr6feCEwpEMEjB5EpJ2lSBAiUFhA+Spff4osJCCDFCm65tQWEkNr1t3oCMwq4cjVjVcyJwLkCrmCd66t3AtMJ5BASbQ8BAgSuEnDqcZW8cQlcK+AE5Fp/oxO4RCCCRw4il0zCoAQIlBYQPkqX3+KLCwggxTeA5dcWEEJq19/qCVwh4MrVFerGJDCXgCtYc9XDbAgMF8ghJNoeAgQInCXg1OMsWf0SWEvACcha9TJbAqcIRPDIQeSUQXRKgEBpAeGjdPktnsAfAgLIHxw+EKgtIITUrr/VEzhDwJWrM1T1SWBtAVew1q6f2RPoLpBDSLQ9BAgQeFfAqce7cn5HYG8BJyB719fqCLwlEMEjB5G3OvEjAgRKCwgfpctv8QS+FBBAvuTxRwK1BYSQ2vW3egLvCLhy9Y6a3xCoJeAKVq16Wy2BlwVyCIm2hwABAs8EnHo8k/E9AQJZwAlI1tAmQOChQASPHEQevuRLAgRKCwgfpctv8QReEhBAXuLyMoHaAkJI7fpbPYFHAq5cPVLxHQECXwm4gvWVjr8RIPBJIIeQaHsIEKgr4NSjbu2tnMBPBJyA/ETPbwkUFYjgkYNIUQbLJlBaQPgoXX6LJ/AjAQHkR3x+TKC2gBBSu/5WX1PAlauadbdqAj0FXMHqqakvAgUFcgiJtocAgX0FnHrsW1srIzBSwAnISG1jEdhUIIJHDiKbLtOyCJQWED5Kl9/iCXQVEEC6cuqMQG0BIaR2/a1+TwFXrvasq1URuFLAFawr9Y1NYEOBHEKi7SFAYF0Bpx7r1s7MCcws4ARk5uqYG4FFBSJ45CCy6DJMm0BpAeGjdPktnsCpAgLIqbw6J1BbQAipXX+rX1PAlas162bWBFYScAVrpWqZK4EFBXIIibaHAIF5BZx6zFsbMyOwk4ATkJ2qaS0EJhWI4JGDyKTTNC0CpQWEj9Llt3gCQwUEkKHcBiNQW0AIqV1/q59TwJWrOetiVgR2FnAFa+fqWhuBCQVyCIm2hwCB6wScelxnb2QClQWcgFSuvrUTuEgggkcOIhdNw7AESgsIH6XLb/EELhUQQC7lNziB2gJCSO36W/01Aq5cXeNuVAIEfgu4gvXbQosAgQsEcgiJtocAgfMEnHqcZ6tnAgTaBZyAtFt5kwCBkwQieOQgctIwuiVQWkD4KF1+iycwlYAAMlU5TIZAbQEhpHb9rf4cAVeuznHVKwEC7wu4gvW+nV8SIHCCQA4h0fYQIPC+gFOP9+38kgCB8wScgJxnq2cCBN4UiOCRg8ib3fgZgdICwkfp8ls8gakFBJCpy2NyBGoLCCG162/17wm4cvWem18RIDBOwBWscdZGIkDgDYEcQqLtIUDguYBTj+c2/kKAwDwCTkDmqYWZECDwRCCCRw4iT17zNYHSAsJH6fJbPIGlBASQpcplsgRqCwghtetv9Y8FXLl67OJbAgTmFXAFa97amBkBAg8EcgiJtodAZQGnHpWrb+0E1hVwArJu7cycQFmBCB45iJSFsPDSAsJH6fJbPIGlBQSQpctn8gRqCwghtetfdfWuXFWtvHUT2EfAFax9amklBEoK5BASbQ+BnQWceuxcXWsjUEfACUidWlspgW0FInjkILLtQi2stIDwUbr8Fk9gKwEBZKtyWgyB2gJCSO3677p6V652rax1Eagr4ApW3dpbOYEtBXIIibaHwMoCTj1Wrp65EyDwTMAJyDMZ3xMgsKxABI8cRJZdiImXFhA+Spff4glsLSCAbF1eiyNQW0AIqV3/VVfvytWqlTNvAgRaBVzBapXyHgECSwrkEBJtD4GZBZx6zFwdcyNAoJeAE5BekvohQGBagQgeOYhMO1ETKy0gfJQuv8UTKCUggJQqt8USqC0ghNSu/6yrd+Vq1sqYFwECZwm4gnWWrH4JEJhSIIeQaHsIXCng1ONKfWMTIHCVgBOQq+SNS4DAZQIRPHIQuWwiBi4tIHyULr/FEygtIICULr/FE6gtIITUrv9Vq3fl6ip54xIgMIuAK1izVMI8CBC4RCCHkGh7CJwp4NTjTF19EyCwioATkFUqZZ4ECJwmEMEjB5HTBtJxaQHho3T5LZ4AgSQggCQMTQIEagsIIbXrf9bqXbk6S1a/BAisKuAK1qqVM28CBE4RyCEk2h4CPxFw6vETPb8lQGBXAScgu1bWuggQeFsggkcOIm935IelBYSP0uW3eAIEvhAQQL7A8ScCBGoLCCG16//u6l25elfO7wgQqCLgClaVSlsnAQJvCeQQEm0Pga8EnHp8peNvBAgQ+LeAExA7gQABAt8IRPDIQeSb1/25qIDwUbTwlk2AwMsCAsjLZH5AgEBVASGkauW/XrcrV1/7+CsBAgSOAq5gHUV8JkCAwBcCOYRE21NbwKlH7fpbPQEC7wk4AXnPza8IECgsEMEjB5HCFKWXLnyULr/FEyDwAwEB5Ad4fkqAQG0BIaRm/V25qll3qyZAoJ+AK1j9LPVEgEBBgRxCou3ZW8Cpx971tToCBMYIOAEZ42wUAgQ2FojgkYPIxkstvTTho3T5LZ4AgY4CAkhHTF0RIFBbQAjZs/6uXO1ZV6siQOA6AVewrrM3MgECGwrkEBJtz9oCTj3Wrp/ZEyAwp4ATkDnrYlYECCwsEMEjB5GFl1J66sJH6fJbPAECJwoIICfi6poAgdoCQsia9Xflas26mTUBAusIuIK1Tq3MlACBBQVyCIm2Z24Bpx5z18fsCBDYQ8AJyB51tAoCBCYWiOCRg8jEUy09NeGjdPktngCBgQICyEBsQxEgUFtACJmz/q5czVkXsyJAYF8BV7D2ra2VESAwoUAOIdH2XCvg1ONaf6MTIFBTwAlIzbpbNQECFwpE8MhB5MKplB5a+ChdfosnQOBCAQHkQnxDEyBQW0AIuab+rlxd425UAgQI3AVcwbpL+CcBAgQuEMghJNqecwWcepzrq3cCBAi0CDgBaVHyDgECBE4UiOCRg8iJQ5XuWvgoXX6LJ0BgIgEBZKJimAoBArUFhJBz6u/K1TmueiVAgMC7Aq5gvSvndwQIEDhBIIeQaHt+JuDU42d+fk2AAIEzBJyAnKGqTwIECPxAIIJHDiI/6Kr0T4WP0uW3eAIEJhYQQCYujqkRIFBbQAh5r/6uXL3n5lcECBAYJeAK1ihp4xAgQOANgRxCou35WsCpx9c+/kqAAIEZBJyAzFAFcyBAgMAXAhE8chD54tXSfxI+Spff4gkQWEhAAFmoWKZKgEBtASHkcf1duXrs4lsCBAjMKuAK1qyVMS8CBAg8EMghJNrVH6ce1XeA9RMgsKKAE5AVq2bOBAiUFojgkYNIVQzho2rlrZsAgdUFBJDVK2j+BAiUFagaQly5KrvlLZwAgU0EXMHapJCWQYBATYEcQqK9++PUY/cKWx8BAhUEnIBUqLI1EiCwtUAEjxxEdl2s8LFrZa2LAIFqAgJItYpbLwEC2wrsGkJcudp2y1oYAQJFBVzBKlp4yyZAYE+BHEKivfrj1GP1Cpo/AQIEPgs4Afls4hsCBAgsLRDBIweRVRcjfKxaOfMmQIDA1wICyNc+/kqAAIFlBVYNIa5cLbvlTJwAAQJNAq5gNTF5iQABAmsK5BAS7dkfpx6zV8j8CBAg8HMBJyA/N9QDAQIEphaI4JGDyKyTFT5mrYx5ESBAoK+AANLXU28ECBCYVmDWEOLK1bRbxsQIECBwioArWKew6pQAAQJzCuQQEu2rH6ceV1fA+AQIEBgv4ARkvLkRCRAgcKlABI8cRK6ajPBxlbxxCRAgcK2AAHKtv9EJECBwmcBVIcSVq8tKbmACBAhMIeAK1hRlMAkCBAhcI5BDSLTPfpx6nC2sfwIECMwv4ARk/hqZIQECBE4ViOCRg8hZgwkfZ8nqlwABAmsJCCBr1ctsCRAgcJrAWSHElavTSqZjAgQILCngCtaSZTNpAgQInCOQQ0i0f/o49fipoN8TIEBgPwEnIPvV1IoIECDwI4EIHjmIvNuZ8PGunN8RIEBgbwEBZO/6Wh0BAgTeFng3hLhy9Ta5HxIgQKCEgCtYJcpskQQIEHhPIIeQaH/3OPX4TsjfCRAgQMAJiD1AgAABAl8KRPDIQeTZy8LHMxnfEyBAgEAWEECyhjYBAgQIPBV4FkJcuXpK5g8ECBAg8EDAFawHKL4iQIAAgccCxxDi1OOxk28JECBA4LmAAPLcxl8I3ATif8HyECDwXMC/Rp7b+AsBAgQIfBYQQD6b+IbAHwLxn/h6CBD4LXAMHP418ttGi8Bd4Pivk/v3/kmAwK9f/jsgdgEBAgQINAvE/1IVgeMeOuKf8Z3/ZauZ0IsECBAoLyCAlN8CAAgQINAmcA8fx7fvgUQIOcr4TIAAAQKPBASQRyq+I0CAAIEPgfsJx/3U4+MPh8b9NOTwtY8ECBAgQOAPAf8dkD84fCBAgACBLPDs1CO/k9s5hHwXWPLvtAkQIECgjoATkDq1tlICBAi8JPBq+Lh3HsEjB5H79/5JgAABAgRCQACxDwgQIEDgD4HWK1d//OjBByHkAYqvCBAgQOCXK1g2AQECBAh8CLx76vHRwaGRQ0i0PQQIECBAwAmIPUCAAAECN4He4ePOGsEjB5H79/5JgAABAjUFBJCadbdqAgQIfAj0unL10eGThhDyBMbXBAgQKCbgClaxglsuAQIEssBZpx55jNzOISTaHgIECBCoJ+AEpF7NrZgAAQI3gdHh484ewSMHkfv3/kmAAAECNQQEkBp1tkoCBAh8CIy6cvUx4JOGEPIExtcECBDYXMAVrM0LbHkECBDIAledeuQ55HYOIdH2ECBAgMD+Ak5A9q+xFRIgQOAmMFv4uJclgkcOIvfv/ZMAAQIE9hQQQPasq1URIEDgQ2CWK1cfE3rSEEKewPiaAAECmwm4grVZQS2HAAECWWDWU488x9zOISTaHgIECBDYT8AJyH41tSICBAjcBFYLH/eyRfDIQeT+vX8SIECAwB4CAsgedbQKAgQIfAiscuXqY8JPGkLIExhfEyBAYHEBV7AWL6DpEyBAIAuseuqR15DbOYRE20OAAAEC6ws4AVm/hlZAgACBm8Bu4eNe1ggeOYjcv/dPAgQIEFhTQABZs25mTYAAgQ+BXa5cfSzoSUMIeQLjawIECCwm4ArWYgUzXQIECGSBXU898hpzO4eQaHsIECBAYD0BJyDr1cyMCRAgcBOoFj7uZY/gkYPI/Xv/JECAAIE1BASQNepklgQIEPgQqHLl6mPBTxpCyBMYXxMgQGByAVewJi+Q6REgQCALVD31yAa5nUNItD0ECBAgML+AE5D5a2SGBAgQuAkIH483QgSPHEQev+VbAgQIEJhFQACZpRLmQYAAgScCrlw9gTl8LYQcQHwkQIDApAKuYE1aGNMiQIBACDj1eG0f5BASbQ8BAgQIzCfgBGS+mpgRAQIEbgLCx3sbIYJHDiLv9eJXBAgQIHCWgABylqx+CRAg8KaAK1dvwh1+JoQcQHwkQIDAJAKuYE1SCNMgQIBACDj16LsPcgiJtocAAQIErhdwAnJ9DcyAAAECNwHh45yNEMEjB5FzRtErAQIECLQKCCCtUt4jQIDASQKuXJ0Ee+hWCDmA+EiAAIGLBFzBugjesAQIEAgBpx5j90EOIdH2ECBAgMB4AScg482NSIAAgZuA8HHNRojgkYPINbMwKgECBOoKCCB1a2/lBAhcJODK1UXwh2GFkAOIjwQIEBgk4ArWIGjDECBAIAScesy1D3IIibaHAAECBM4XcAJyvrERCBAgcBMQPubcCBE8chCZc5ZmRYAAgX0EBJB9amklBAhMKuDK1aSFOUxLCDmA+EiAAIGTBFzBOglWtwQIEAgBpx5r7YMcQqLtIUCAAIH+Ak5A+pvqkQABAjcB4WPNjRDBIweRNVdh1gQIEJhXQACZtzZmRoDAogKuXC1auMO0hZADiI8ECBDoJOAKVidI3RAgQCAEnHrstQ9yCIm2hwABAgR+LuAE5OeGeiBAgMBNQPjYcyNE8MhBZM9VWhUBAgTGCQgg46yNRIDApgKuXG1a2MOyhJADiI8ECBB4U8AVrDfh/IwAAQIh4NSj1j7IISTaHgIECBB4XcAJyOtmfkGAAIGbgPBRcyNE8MhBpKaCVRMgQOB9AQHkfTu/JECgqIArV0ULf1i2EHIA8ZEAAQKNAq5gNUJ5jQABAiHg1MM+yAI5hETbQ4AAAQLfCzgB+d7IGwQIELgJCB82wiOBCB45iDx6x3cECBAg8FtAAPltoUWAAIGHAq5cPWTx5UFACDmA+EiAAIEnAq5gPYHxNQECBELAqYd98IpADiHR9hAgQIDAZwEnIJ9NfEOAAIGbgPBhI7wjEMEjB5F3+vAbAgQI7CwggOxcXWsjQOAtAVeu3mLzo4OAEHIA8ZEAAQJ/CbiCZSsQIEAgCTj1SBiaPxbIISTaHgIECBD49csJiF1AgACBvwSED1vhDIEIHjmInDGGPgkQILCSgACyUrXMlQCBUwRcuTqFVacHASHkAOIjAQJlBVzBKlt6CydAIAScetgHIwVyCIm2hwABAhUFnIBUrLo1EyBwExA+bIQrBCJ45CByxRyMSYAAgSsFBJAr9Y1NgMAlAq5cXcJu0IOAEHIA8ZEAgTICrmCVKbWFEiAQAk497IOZBHIIibaHAAECFQScgFSosjUSIHATED5shBkFInjkIDLjHM2JAAECPQUEkJ6a+iJAYEoBV66mLItJHQSEkAOIjwQIbCvgCta2pbUwAgRCwKmHfbCSQA4h0fYQIEBgRwEnIDtW1ZoIELgJCB82wooCETxyEFlxDeZMgACBrwQEkK90/I0AgSUFXLlasmwmfRAQQg4gPhIgsI2AK1jblNJCCBAIAace9sFOAjmERNtDgACBHQScgOxQRWsgQOAmIHzYCDsKRPDIQWTHNVoTAQK1BASQWvW2WgJbCrhytWVZLeogIIQcQHwkQGBZAVewli2diRMgEAJOPeyDSgI5hETbQ4AAgRUFnICsWDVzJkDgJiB82AgVBSJ45CBS0cCaCRBYW0AAWbt+Zk+gpIArVyXLbtEHASHkAOIjAQLLCLiCtUypTJQAgRBw6mEfEPgtkENItD0ECBBYQcAJyApVMkcCBG4CwoeNQOCzQASPHEQ+v+EbAgQIzCUggMxVD7MhQOCBgCtXD1B8ReAgIIQcQHwkQGBaAVewpi2NiREgEAJOPewDAu0COYRE20OAAIEZBZyAzFgVcyJA4CYgfNgIBF4XiOCRg8jrPfgFAQIEzhUQQM711TsBAm8IuHL1BpqfEDgICCEHEB8JEJhGwBWsaUphIgQIhIBTD/uAQD+BHEKi7SFAgMAMAk5AZqiCORAgcBMQPmwEAv0FInjkINJ/BD0SIEDgNQEB5DUvbxMgcIKAK1cnoOqSwEFACDmA+EiAwGUCrmBdRm9gAgRCwKmHfUBgnEAOIdH2ECBA4AoBJyBXqBuTAIGbgPBhIxAYLxDBIweR8TMwIgEC1QUEkOo7wPoJXCDgytUF6IYkcBAQQg4gPhIgMEzAFaxh1AYiQCAEnHrYBwTmEcghJNoeAgQIjBBwAjJC2RgECNwEhA8bgcB8AhE8chCZb4ZmRIDAbgICyG4VtR4CEwq4cjVhUUyJwEFACDmA+EiAwGkCrmCdRqtjAgRCwKmHfUBgHYEcQqLtIUCAwBkCTkDOUNUnAQI3AeHDRiCwnkAEjxxE1luBGRMgMLuAADJ7hcyPwIICrlwtWDRTJnAQEEIOID4SINBNwBWsbpQ6IkAgBJx62AcE9hHIISTaHgIECPQQcALSQ1EfBAjcBIQPG4HAfgIRPHIQ2W+FVkSAwGgBAWS0uPEIbCjgytWGRbUkAgcBIeQA4iMBAm8LuIL1Np0fEiAQAk497AMCdQRyCIm2hwABAu8IOAF5R81vCBC4CQgfNgKBegIRPHIQqSdgxQQI/FRAAPmpoN8TKCjgylXBolsygYOAEHIA8ZEAgWYBV7CaqbxIgEAIOPWwDwgQuAvkEBJtDwECBFoEnIC0KHmHAIGbgPBhIxAgcBSI4JGDyPHvPhMgQOAoIIAcRXwmQOCTgCtXn0h8QYDAQUAIOYD4SIDAUwFXsJ7S+AMBAiHg1MM+IECgVSCHkGh7CBAg8EjACcgjFd8RIHATED5sBAIEXhWI4JGDyKu/9z4BAvsLCCD719gKCbws4MrVy2R+QIDAQUAIOYD4SIDAh4ArWB8UGgQIhIBTD/uAAIFeAjmERNtDgACBEHACYh8QIPAhIHx8UGgQINBJIIJHDiKdutUNAQILCwggCxfP1An0EnDlqpekfggQeCYghDyT8T2BegKuYNWruRUT+EPAqccfHD4QIHCiQA4h0fYQIFBTwAlIzbpbNYGbgPBhIxAgMFoggkcOIqPHNx4BAtcLCCDX18AMCAwXcOVqOLkBCRA4CAghBxAfCRQScAWrULEtlUAIOPWwDwgQmEUgh5BoewgQqCHgBKRGna2SwE1A+LARCBCYTSCCRw4is83PfAgQ6C8ggPQ31SOB6QRcuZquJCZEgMBBQAg5gPhIYGMBV7A2Lq6lEQgBpx72AQECqwjkEBJtDwECewo4AdmzrlZF4CYgfNgIBAisJhDBIweR1eZvvgQIfC8ggHxv5A0Cywm4crVcyUyYAIGDgBByAPGRwEYCrmBtVExLIRACTj3sAwIEdhHIISTaHgIE9hBwArJHHa2CwE1A+LARCBDYTSCCRw4iu63PeghUFBBAKlbdmrcTcOVqu5JaEAECBwEh5ADiI4GFBVzBWrh4pk4gBJx62AcECFQRyCEk2h4CBNYUcAKyZt3MmsBNQPiwEQgQqCYQwSMHkWrrt14COwgIIDtU0RrKCbhyVa7kFkyAwEFACDmA+EhgIQFXsBYqlqkSCAGnHvYBAQIE/i2QQ0i0PQQIrCHgBGSNOpklgZuA8GEjECBA4E+BCB45iPz5V58IEJhRQACZsSrmROAg4MrVAcRHAgQIHASEkAOIjwQmFnAFa+LimBqBEHDqYR8QIECgTSCHkGh7CBCYU8AJyJx1MSsCNwHhw0YgQIDAawIRPHIQee3X3iZAYISAADJC2RgEXhRw5epFMK8TIEDgICCEHEB8JDCRgCtYExXDVAiEgFMP+4AAAQJ9BHIIibaHAIE5BJyAzFEHsyBwExA+bAQCBAj0FYjgkYNI3971RoDAOwICyDtqfkOgs4ArV51BdUeAAIGDgBByAPGRwIUCrmBdiG9oAiHg1MM+IECAwBiBHEKi7SFA4BoBJyDXuBuVwE1A+LARCBAgMFYggkcOImNHNxoBAiEggNgHBC4QcOXqAnRDEiBAIAkIIQlDk8BgAVewBoMbjoBTD3uAAAECcwjkEBJtDwECYwScgIxxNgqBm4DwYSMQIEBgLoEIHjmIzDU7syGwp4AAsmddrWoyAVeuJiuI6RAgQOAgIIQcQHwkcKKAK1gn4uqaQAg49bAPCBAgsIZADiHR9hAgcI6AE5BzXPVK4CYgfNgIBAgQWEsggkcOImvN3mwJrCEggKxRJ7NcTMCVq8UKZroECBA4CAghBxAfCXQUcAWrI6auCISAUw/7gAABAnsI5BASbQ8BAn0EnID0cdQLgZuA8GEjECBAYC+BCB45iOy1OqshcI2AAHKNu1E3E3DlarOCWg4BAgQOAkLIAcRHAj8QcAXrB3h+SiAEnHrYBwQIEKghkENItD0ECLwn4ATkPTe/InATED5sBAIECNQSiOCRg0it1VstgT4CTkD6OOqlmEAEj3j8J2DFCm+5BAgQ+EvgHkI+/fvAP//9wt9//euE5L//evlv2AgQyAICSNbQJtAg4NSjAckrBAgQKCBwDyGx1FsQ+St8fFp6fC+EfGLxRV0BAaRu7a38DQHh4w00PyFAgMDGAh8nIM/Cx33tQshdwj8J/BJAbAICDQKuXDUgeYUAAQJVBb4LH3cXIeQu4Z/FBQSQ4hvA8r8XcOrxvZE3CBAgQIAAAQKtAnEjsTW3t/bpPQLbCNxPPrZZkIUQIECAQHeBv//3C/8nef13Qbr763A9AQFkvZqZMQECBAgQIDCTwCv/Ua4AMlPlzOUiAf9/QC6CNywBAgQIECBAgACBigICSMWqWzMBAgQIECDQT6D1VKP1vX4z0xOBKQUEkCnLYlIECBAgQIDAUgLfhYvv/r7UYk2WwM8EBJCf+fk1AQIECBAgQODfAs9CxrPvuREoKhD/knjlvzpVlMmyCRAgQIAAAQIECBDoIeAEpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBADwEBpIeiPggQIECAAAECBAgQaBIQQJqYvESAAAECBAgQIECAQA8BAaSHoj4IECBAgAABAgQIEGgSEECamLxEgAABAgQIECBAgEAPAQGkh6I+CBAgQIAAAQIECBBoEhBAmpi8RIAAAQIECBAgQIBAD4H/B3GXoHtV+iqwAAAAAElFTkSuQmCC", - "text/html": [ - "" - ], - "text/plain": [ - " Size: 2MB\n", - "array([[4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " ...,\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080]], shape=(600, 800), dtype=uint32)\n", - "Coordinates:\n", - " * y (y) float64 5kB 0.003333 0.01 0.01667 0.02333 ... 3.983 3.99 3.997\n", - " * x (x) float64 6kB 0.0025 0.0075 0.0125 0.0175 ... 3.987 3.993 3.998" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "from xrspatial import a_star_search\n", - "\n", "# define range of x and y\n", "xrange = (0, 4)\n", "yrange = (0, 4)\n", "\n", "# create line raster\n", "ys = [1, 1, 3, 3, 1, 1, np.nan, 1, 3, np.nan, 1, 3, np.nan, 1, 3, np.nan, 2, 2]\n", - "xs = [1, 3, 3, 1, 1, 3, np.nan, 1, 3, np.nan, 3, 1, np.nan, 2, 2, np.nan, 1, 3]\n", + "xs = [1, 3, 3, 1, 1, 3, np.nan, 1, 3, np.nan, 3, 1, np.nan, 1, 2, np.nan, 1, 3]\n", "line_df = pd.DataFrame(dict(x=xs, y=ys))\n", "\n", "W = 800\n", @@ -143,7 +60,7 @@ "line_agg = cvs.line(line_df, x=\"x\", y=\"y\").astype(int)\n", "line_shaded = dynspread(shade(line_agg, cmap=[\"black\", \"salmon\"]))\n", "\n", - "# pick up 2 random locations\n", + "# pick start and goal locations\n", "start = (3, 1)\n", "goal = (1, 3)\n", "\n", @@ -162,80 +79,44 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We're now ready to apply `a_star_search`." + "#### 8-connectivity shortest path\n", + "\n", + "Cells with value 0 (off the lines) are barriers. The path is highlighted in green." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "#### Calculating the 8-connectivity shortest path" + "path_agg_8 = a_star_search(\n", + " line_agg, start, goal, barriers=[0], snap_start=True, snap_goal=True\n", + ")\n", + "\n", + "path_shaded = dynspread(shade(path_agg_8, cmap=[\"green\"]))\n", + "set_background(stack(line_shaded, path_shaded, start_shaded, goal_shaded), \"black\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "- To calculate the path, we input the line raster and the start and goal point coordinates.\n", - "- We also set the barriers; i.e. cells that are not crossable. In our case, any cell with a value of 0 (all the black non-line cells).\n", - "- Finally, we'll also set snap_start and snap_goal to True.\n", - "- Note: since `a_star_search` uses 8-connectivity by default, we don't need to pass that in.\n", - "\n", - "The shortest path is highlighted in the rendering below." + "#### 4-connectivity shortest path" ] }, { "cell_type": "code", - "execution_count": 3, - "metadata": { - "execution": { - "iopub.execute_input": "2026-01-23T01:07:49.995767Z", - "iopub.status.busy": "2026-01-23T01:07:49.995587Z", - "iopub.status.idle": "2026-01-23T01:07:50.933357Z", - "shell.execute_reply": "2026-01-23T01:07:50.932684Z" - } - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAA4fklEQVR4Ae3da5LkLJIo0OywsW/Bs4HxFfSCZ370gCqoqCwpMx4ZAgk/MoOKRA/wA23XvKHn/uvj4+M/pbgIECBAgAABAgQIECCwu8Bl9x50QIAAAQIECBAgQIAAgauABMRSIECAAAECBAgQIECgm4AEpBu1jggQIECAAAECBAgQkIBYAwQIECBAgAABAgQIdBOQgHSj1hEBAgQIECBAgAABAhIQa4AAAQIECBAgQIAAgW4CEpBu1DoiQIAAAQIECBAgQEACYg0QIECAAAECBAgQINBNQALSjVpHBAgQIECAAAECBAhIQKwBAgQIECBAgAABAgS6CUhAulHriAABAgQIECBAgAABCYg1QIAAAQIECBAgQIBANwEJSDdqHREgQIAAAQIECBAgIAGxBggQIECAAAECBAgQ6CYgAelGrSMCBAgQIECAAAECBCQg1gABAgQIECBAgAABAt0EJCDdqHVEgAABAgQIECBAgIAExBogQIAAAQIECBAgQKCbgASkG7WOCBAgQIAAAQIECBCQgFgDBAgQIECAAAECBAh0E5CAdKPWEQECBAgQIECAAAECEhBrgAABAgQIECBAgACBbgISkG7UOiJAgAABAgQIECBAQAJiDRAgQIAAAQIECBAg0E1AAtKNWkcECBAgQIAAAQIECEhArAECBAgQIECAAAECBLoJSEC6UeuIAAECBAgQIECAAAEJiDVAgAABAgQIECBAgEA3AQlIN2odESBAgAABAgQIECAgAbEGCBAgQIAAAQIECBDoJiAB6UatIwIECBAgQIAAAQIEJCDWAAECBAgQIECAAAEC3QQkIN2odUSAAAECBAgQIECAgATEGiBAgAABAgQIECBAoJuABKQbtY4IECBAgAABAgQIEJCAWAMECBAgQIAAAQIECHQTkIB0o9YRAQIECBAgQIAAAQISEGuAAAECBAgQIECAAIFuAhKQbtQ6IkCAAAECBAgQIEBAAmINECBAgAABAgQIECDQTUAC0o1aRwQIECBAgAABAgQISECsAQIECBAgQIAAAQIEuglIQLpR64gAAQIECBAgQIAAAQmINUCAAAECBAgQIECAQDcBCUg3ah0RIECAAAECBAgQICABsQYIECBAgAABAgQIEOgmIAHpRq0jAgQIECBAgAABAgQkINYAAQIECBAgQIAAAQLdBCQg3ah1RIAAAQIECBAgQICABMQaIECAAAECBAgQIECgm4AEpBu1jggQIECAAAECBAgQkIBYAwQIECBAgAABAgQIdBOQgHSj1hEBAgQIECBAgAABAhIQa4AAAQIECBAgQIAAgW4CEpBu1DoiQIAAAQIECBAgQEACYg0QIECAAAECBAgQINBNQALSjVpHBAgQIECAAAECBAhIQKwBAgQIECBAgAABAgS6CUhAulHriAABAgQIECBAgAABCYg1QIAAAQIECBAgQIBANwEJSDdqHREgQIAAAQIECBAgIAGxBggQIECAAAECBAgQ6CYgAelGrSMCBAgQIECAAAECBCQg1gABAgQIECBAgAABAt0EJCDdqHVEgAABAgQIECBAgIAExBogQIAAAQIECBAgQKCbgASkG7WOCBAgQIAAAQIECBCQgFgDBAgQIECAAAECBAh0E5CAdKPWEQECBAgQIECAAAECEhBrgAABAgQIECBAgACBbgISkG7UOiJAgAABAgQIECBAQAJiDRAgQIAAAQIECBAg0E1AAtKNWkcECBAgQIAAAQIECEhArAECBAgQIECAAAECBLoJSEC6UeuIAAECBAgQIECAAAEJiDVAgAABAgQIECBAgEA3AQlIN2odESBAgAABAgQIECAgAbEGCBAgQIAAAQIECBDoJiAB6UatIwIECBAgQIAAAQIEJCDWAAECBAgQIECAAAEC3QQkIN2odUSAAAECBAgQIECAgATEGiBAgAABAgQIECBAoJuABKQbtY4IECBAgAABAgQIEJCAWAMECBAgQIAAAQIECHQTkIB0o9YRAQIECBAgQIAAAQISEGuAAAECBAgQIECAAIFuAhKQbtQ6IkCAAAECBAgQIEBAAmINECBAgAABAgQIECDQTUAC0o1aRwQIECBAgAABAgQISECsAQIECBAgQIAAAQIEuglIQLpR64gAAQIECBAgQIAAAQmINUCAAAECBAgQIECAQDcBCUg3ah0RIECAAAECBAgQICABsQYIECBAgAABAgQIEOgmIAHpRq0jAgQIECBAgAABAgQkINYAAQIECBAgQIAAAQLdBCQg3ah1RIAAAQIECBAgQICABMQaIECAAAECBAgQIECgm4AEpBu1jggQIECAAAECBAgQkIBYAwQIECBAgAABAgQIdBOQgHSj1hEBAgQIECBAgAABAhIQa4AAAQIECBAgQIAAgW4CEpBu1DoiQIAAAQIECBAgQEACYg0QIECAAAECBAgQINBNQALSjVpHBAgQIECAAAECBAhIQKwBAgQIECBAgAABAgS6CUhAulHriAABAgQIECBAgAABCYg1QIAAAQIECBAgQIBANwEJSDdqHREgQIAAAQIECBAgIAGxBggQIECAAAECBAgQ6CYgAelGrSMCBAgQIECAAAECBCQg1gABAgQIECBAgAABAt0EJCDdqHVEgAABAgQIECBAgIAExBogQIAAAQIECBAgQKCbgASkG7WOCBAgQIAAAQIECBCQgFgDBAgQIECAAAECBAh0E5CAdKPWEQECBAgQIECAAAECEhBrgAABAgQIECBAgACBbgISkG7UOiJAgAABAgQIECBAQAJiDRAgQIAAAQIECBAg0E1AAtKNWkcECBAgQIAAAQIECEhArAECBAgQIECAAAECBLoJSEC6UeuIAAECBAgQIECAAAEJiDVAgAABAgQIECBAgEA3AQlIN2odESBAgAABAgQIECAgAbEGCBAgQIAAAQIECBDoJiAB6UatIwIECBAgQIAAAQIEJCDWAAECBAgQIECAAAEC3QQkIN2odUSAAAECBAgQIECAgATEGiBAgAABAgQIECBAoJuABKQbtY4IECBAgAABAgQIEJCAWAMECBAgQIAAAQIECHQTkIB0o9YRAQIECBAgQIAAAQISEGuAAAECBAgQIECAAIFuAhKQbtQ6IkCAAAECBAgQIEBAAmINECBAgAABAgQIECDQTUAC0o1aRwQIECBAgAABAgQISECsAQIECBAgQIAAAQIEuglIQLpR64gAAQIECBAgQIAAAQmINUCAAAECBAgQIECAQDcBCUg3ah0RIECAAAECBAgQICABsQYIECBAgAABAgQIEOgmIAHpRq0jAgQIECBAgAABAgT+CwEBAmuB/6ybPv610aaJAAECBAgQIEDgOQE7IM95eTqBwFbyUcP+qj0BiRAJECBAgAABAm8TkIC8jdKHZhC4l2Tcuz+DgRgIECBAgAABAnsKSED21PXtUwk8mlw8+typgjdYAgQIECBAgEAnAQlIJ2jdECBAgAABAgQIECDw8eF/hG4VEHhB4H/jv194yysECBAgkEngn/h3pnDFSuBhAQnIw1QeJHAT8P+o3Cz8yitQE3H/Wcg7/yLfFvBfUG27aCXwp4AjWH9q+E2AAAECBAgQeFGgJeUS8xcBvZZGQAKSZqoFek/g0f9/Ph597l5/7hMgQIDAPAIt+ZgnIpEQ2E/AEaz9bH35hAI1ufju/8rVknxEeagWFwECBAikF2hHrux6pF8KAJ4QkIA8geXRHAJfJSG/dz6iONRSr1hqFQECBAgkFLDrkXDShfwWAUew3sLoI7MJ1GSjlv8r/yPb9vtTjFH+auXTDX8QIECAQAYByUeGWRbjXgJ2QPaS9d0cAlHCbCVHxKIkQIBAagFHrlJPv+DfJCABeROkzyQWiBJ7LfWKpVYRIECAwIQCdj0mnFQhDRG4DOlVpwRmE4gSUCuzxSYeAgQIEPiQfFgEBN4nYAfkfZa+ROCWhAQMAgQIEJhBwJGrGWZRDEcTkIAcbUaM5/wCUUKopV6x1CoCBAgQOKGAXY8TTpohn0LAEaxTTJNBnk4gyohbOd3gDZgAAQIEJB/WAIH9BOyA7GfrywRuSUjAIECAAIEzCDhydYZZMsazC0hAzj6Dxn98gShDrKVesdQqAgQIEDiggF2PA06KIU0p4AjWlNMqqMMJRBlRK4cbnAERIECAgOTDGiDQT8AOSD9rPRG4JSEBgwABAgSOIODI1RFmwRiyCUhAss24eMcLRBlCLfWKpVYRIECAwAABux4D0HVJoAg4gmUZEBghEKXTVkb0r08CBAgkF5B8JF8Awh8qYAdkKL/O0wtEEWglPQYAAgQI7C/gyNX+xnogcE9AAnJPyH0CewtE6aCWesVSqwgQIEBgBwG7Hjug+iSBFwQuL7zjFQIE3i0Q5YOtvPvbvkeAAAECH5IPi4DAcQTsgBxnLoyEwC0JCRgECBAg8A4BR67eoegbBN4rIAF5r6evEfi5QJRP1FKvWGoVAQIECLwgYNfjBTSvEOgg4AhWB2RdEHhaIMobrTz9shcIECBAQPJhDRA4roAdkOPOjZERuCUhAYMAAQIEHhFw5OoRJc8QGCsgARnrr3cC9wWiPFJLvWKpVQQIECCwIWDXYwNFE4EDCjiCdcBJMSQCK4EoLa2sbmogQIAAAcmHNUDgPAJ2QM4zV0ZK4JaEBAwCBAgQqAKOXFkHBM4nIAE535wZcXaBKAC11CuWWkWAAIGUAnY9Uk67oCcQcARrgkkUQkKBKDG3kjB8IRMgQEDyYQ0QOK+AHZDzzp2RE7glIQGDAAECOQQcucoxz6KcW0ACMvf8ii6DQJQga6lXLLWKAAECUwrY9ZhyWgWVUMARrISTLuQJBaLE1MqE4QmJAAECkg9rgMA8AnZA5plLkRC4JSEBgwABAnMIOHI1xzyKgsCfAhKQPzX8JjCDQJQgaqlXLLWKAAECpxSw63HKaTNoAncFHMG6S+QBAicUiDLmVk44fEMmQICA5MMaIDCvgB2QeedWZARuSUjAIECAwDkEHLk6xzwZJYGfCEhAfqLnXQJnEIgyyFrqFUutIkCAwCEF7HoccloMisDbBRzBejupDxI4oECUMbVywOEZEgECBCQf1gCBPAJ2QPLMtUgJ3JKQgEGAAIFjCDhydYx5MAoCPQUkID219UXgCAJRBlFLvWKpVQQIEBgiYNdjCLtOCQwXcARr+BQYAIEBAlH6bGVA97okQICA5MMaIJBXwA5I3rkXOYFbEhIwCBAg0EfAkas+znohcGQBCciRZ8fYCPQQiNJJLfWKpVYRIEBgFwG7Hruw+iiB0wk4gnW6KTNgAjsIRPlmKzt83icJECAg+bAGCBBoAnZAmoR/CRC4JSEBgwABAu8RcOTqPY6+QmAmAQnITLMpFgLvEIjykVrqFUutIkCAwEsCdj1eYvMSgekFHMGafooFSOAFgSjvtPLC614hQICA5MMaIEDgKwE7IF/JaCdA4JaEBAwCBAg8JuDI1WNOniKQWUACknn2xU7gEYEoD9VSr1hqFQECBDYF7HpssmgkQOAvAUew/gLxJwECGwJR2lrZuK2JAAECkg9rgACBRwXsgDwq5TkCBG5JSMAgQIDALwFHrqwEAgSeFZCAPCvmeQLZBaIA1FKvWGoVAQJJBex6JJ14YRP4oYAjWD8E9DqBlAJRom4lJYCgCRCQfFgDBAi8KmAH5FU57xEgcEtCAgYBAlkEHLnKMtPiJLCfgARkP1tfJpBDIEqYtdQrllpFgMCkAnY9Jp1YYRHoLOAIVmdw3RGYUiBKVK1MGaCgCBCQfFgDBAi8S8AOyLskfYcAgVsSEjAIEJhFwJGrWWZSHASOIyABOc5cGAmBOQSihFFLvWKpVQQInFTArsdJJ86wCRxcwBGsg0+Q4RE4pUCUUbdyygAMmgAByYc1QIDAXgJ2QPaS9V0CBG5JSMAgQOAsAo5cnWWmjJPAeQUkIOedOyMncA6BKMOspV6x1CoCBA4qYNfjoBNjWAQmE3AEa7IJFQ6BQwpEGVUrhxygQREgIPmwBggQ6CVgB6SXtH4IELglIQGDAIGjCDhydZSZMA4CeQQkIHnmWqQEjiEQZRi11CuWWkWAwCABux6D4HVLILmAI1jJF4DwCQwRiNJrK0MGoFMCBCQf1gABAqME7ICMktcvAQK3JCRgECDQS8CRq17S+iFA4CsBCchXMtoJEOgjEKWbWuoVS60iQGAnAbseO8H6LAECTwk4gvUUl4cJENhFIMpXW9mlAx8lQEDyYQ0QIHAUATsgR5kJ4yBA4JaEBAwCBN4l4MjVuyR9hwCBdwlIQN4l6TsECLxHIMpnaqlXLLWKAIEXBex6vAjnNQIEdhVwBGtXXh8nQOAlgShvtfLSB7xEgIDkwxogQOCoAnZAjjozxkWAwC0JCRgECDwq4MjVo1KeI0BglIAEZJS8fgkQeEwgymO11CuWWkWAwBcCdj2+gNFMgMChBBzBOtR0GAwBApsCUVpb2XxAIwECkg9rgACBswjYATnLTBknAQK3JCRgECDQBBy5ahL+JUDgLAISkLPMlHESIPBLIMo/tdQrllpFIK2AXY+0Uy9wAqcWcATr1NNn8ASSCkSJu5WkBMImIPmwBggQOKuAHZCzzpxxEyBwS0ICBoE8Ao5c5ZlrkRKYVUACMuvMiotAFoEogdZSr1hqFYFpBex6TDu1AiOQSsARrFTTLVgCkwpEiauVSUMUFgHJhzVAgMAsAnZAZplJcRAgcEtCAgaBeQQcuZpnLkVCgMAvAQmIlUCAwFwCUcKppV6x1CoCpxWw63HaqTNwAgS+EXAE6xsctwgQOKlAlHG3ctIQDJuA5MMaIEBgVgE7ILPOrLgIELglIQGDwHkEHLk6z1wZKQECrwlIQF5z8xYBAmcRiDLQWuoVS60icFgBux6HnRoDI0DgjQKOYL0R06cIEDioQJRxtXLQIRoWAcmHNUCAQBYBOyBZZlqcBAjckpCAQeA4Ao5cHWcujIQAgT4CEpA+znohQOAoAlEGUku9YqlVBIYJ2PUYRq9jAgQGCjiCNRBf1wQIDBKI0m8rg4agWwKSD2uAAIGsAnZAss68uAkQuCUhAYNAPwFHrvpZ64kAgWMKSECOOS9GRYBAL4EoHdVSr1hqFYHdBOx67EbrwwQInEjAEawTTZahEiCwk0CU77ayUxc+S0DyYQ0QIEDgl4AdECuBAAECTSDKj1Zam38J/FDAkasfAnqdAIHpBCQg002pgAgQ+JFAlLdrqVcstYrAywJ2PV6m8yIBAhMLXCaOTWgECBB4TSDKa6289gVvEfiQfFgEBAgQ2BawA7LtopUAAQK3JCRgEHhcwJGrx608SYBATgEJSM55FzUBAo8KRHmwlnrFUqsIfClg1+NLGjcIECDwW8ARrN8UfhAgQOALgSjtrXzxiGYCkg9rgAABAo8J2AF5zMlTBAgQuCUhAYPATcCRq5uFXwQIEHhEQALyiJJnCBAg0ASi/KilXrHUqsQCdj0ST77QCRB4WcARrJfpvEiAQFqBKJG3khZB4JIPa4AAAQKvCdgBec3NWwQIELglIQEjk4AjV5lmW6wECOwhIAHZQ9U3CRDIIxAl1FrqFUutmljArsfEkys0AgS6CTiC1Y1aRwQITCsQJbJWpg1SYJIPa4AAAQLvEbAD8h5HXyFAgMAtCQkYMwk4cjXTbIqFAIEjCEhAjjALxkCAwDwCUUKppV6x1KoTC9j1OPHkGToBAocVcATrsFNjYAQInFYgyshbOW0QBi75sAYIECCwj4AdkH1cfZUAAQK3JCRgnEnAkaszzZaxEiBwRgEJyBlnzZgJEDiPQJSh1lKvWGrVgQXsehx4cgyNAIFpBBzBmmYqBUKAwGEFooyslcMO0sAkH9YAAQIE+gjYAenjrBcCBAjckpCAcSQBR66ONBvGQoBABgEJSIZZFiMBAscRiDKUWuoVS60aKGDXYyC+rgkQSCvgCFbaqRc4AQLDBKL03MqwQehY8mENECBAYIyAHZAx7nolQIDALQkJGD0FHLnqqa0vAgQIrAUkIGsTLQQIEOgnEKWrWuoVS63aUcCux464Pk2AAIEHBRzBehDKYwQIENhNIMqXW9mtEx+WfFgDBAgQOIaAHZBjzINRECBA4JaEBIx3Cjhy9U5N3yJAgMDPBSQgPzf0BQIECLxPIMqnaqlXLLXqBwJ2PX6A51UCBAjsJOAI1k6wPkuAAIGXBaK82crLH/Gi5MMaIECAwDEF7IAcc16MigABArckJGA8I+DI1TNaniVAgEB/AQlIf3M9EiBA4HGBKI/WUq9YatU3AnY9vsFxiwABAgcRcATrIBNhGAQIEPhSIMqdVr58yA3JhzVAgACBcwjYATnHPBklAQIEbklIwPhTwJGrPzX8JkCAwPEFJCDHnyMjJECAwE0gys9a6hVLnbqy65F6+gVPgMBJBRzBOunEGTYBAokFosTeSmIGyUfiyRc6AQKnFrADcurpM3gCBFILRIm+lUQQjlwlmmyhEiAwpYAEZMppFRQBAmkEokRaS71iqaeu7HpMPb2CI0AgicAlSZzCJECAwLwCUUJrZd4oPyQfE0+u0AgQSCVgByTVdAuWAIGpBaJE18pEgTpyNdFkCoUAAQJFQAJiGRAgQGAmgSjB1FKvWOpTV3Y9Tj19Bk+AAIFNgctmq0YCBAgQOK9AlKG3ct4oHLk68dwZOgECBL4TsAPynY57BAgQOLNAlMG3cqI4HLk60WQZKgECBF4QkIC8gOYVAgQInEYgykhrqVcs9aErR64OPT0GR4AAgbcIXN7yFR8hQIAAgeMKRBlaK8cdpSNXB54bQyNAgMA7BeyAvFPTtwgQIHBkgSiDa+VA43Tk6kCTYSgECBDoICAB6YCsCwIECBxGIMpIaqlXLPXQypGrofw6J0CAwBCBy5BedUqAAAEC4wSidN3KuFE4cjXQXtcECBAYKWAHZKS+vgkQIDBSIErnrXQchyNXHbF1RYAAgQMKSEAOOCmGRIAAgW4CUXqqpV6x1LtWjlztyuvjBAgQOIXA5RSjNEgCBAgQ2E8gyqdb2a8XR652tPVpAgQInEnADsiZZstYCRAgsKdAlI+38sZ+HLl6I6ZPESBAYAIBCcgEkygEAgQIvE0gypdqqVcs9Y8qR65+xOdlAgQITClwmTIqQREgQIDA6wJRXm3l9a84cvUDO68SIEBgZgE7IDPPrtgIECDwE4EoL7fyxHccuXoCy6MECBBIKCABSTjpQiZAgMDDAlGerKVesdTfVo5cfcvjJgECBAgUgQsFAgQIECDwrUCUu61886Dk4xsctwgQIEDgt4AdkN8UfhAgQIDAtwJR7rbyx4OOXP2B4ScBAgQI3BWQgNwl8gABAgQI/BaI8quW62XXo0n4lwABAgQeFZCAPCrlubQC7b/dTQsgcAJfCPzz8e+P//347+X/2tUXj2gmQIAAAQIrAQnIikQDgc8C/8S/Pzf4i0BygZaU1+SjJiF/7ogkpxE+gd8C7T8nvxv8IEDgt8Dl9y8/CBAgQIDAHYF25Op3Yh7lhVbuvOs2AQIECBCoAhIQ64AAAQIEHhJoycfq4SgtraxuaiBAgAABAp8FHMH67OEvAgQIEPhLoB0l+b3r8df9339G+dXK70Y/CBAgQIDAZwEJyGcPfxEgQIDAHwJf7nr88cynn1H+qqVesdQqAgQIECDwSeDy6S9/ECBAgACBq8DTyUeTi/KjldbmXwIECBAgcBWwA2IpECBAgMAngYePXH16a+OPKG2tbNzWRIAAAQI5BSQgOedd1AQIENgUeHnXY/NrpTGupd6vv10ECBAgkF7gkl4AAAECBAgsAm9PPpprlB+ttDb/EiBAgEBaATsgaade4AQIEPgl8LYjV/dAozzQyr1n3SdAgACBaQUkINNOrcAIECBwX2C3XY+vuo5yo5Z6xVKrCBAgQCCZwCVZvMIlQIAAgatA9+SjyUf50Upr8y8BAgQIpBGwA5JmqgVKgACBXwLdjlzdA4/yQCv3nnWfAAECBKYRkIBMM5UCIUCAwH2BYbseXw0tyo1a6hVLrSJAgACByQUuk8cnPAIECBC4Chwu+WgzE+VHK63NvwQIECAwrYAdkGmnVmAECBD4JXCYI1f3JiTKA63ce9Z9AgQIEDitgATktFNn4AQIELgvcNhdj6+GHuVGLfWKpVYRIECAwGQCl8niEQ4BAgQIXAVOl3y0mYvyo5XW5l8CBAgQmEbADsg0UykQAgQI/BI4zZGrexMW5YFW7j3rPgECBAicRkACcpqpMlACBAjcFzjtrsdXoUW5UUu9YqlVBAgQIHBygcvJx2/4BAgQIHAVmC75aDMb5Ucrrc2/BAgQIHBaATsgp506AydAgMAvgWmOXN2b0CgPtHLvWfcJECBA4LACEpDDTo2BESBA4L7AtLseX4Ue5UYt9YqlVhEgQIDAyQQuJxuv4RIgQIDAVSBd8tFmPsqPVlqbfwkQIEDgNAJ2QE4zVQZKgACBXwJpjlzdm/AoD7Ry71n3CRAgQOAwAhKQw0yFgRAgQOC+QNpdj69ootyopV6x1CoCBAgQOLjA5eDjMzwCBAgQuApIPr5YClHaW/niEc0ECBAgcBwBOyDHmQsjIUCAwKaAI1ebLOvGKE2trO9qIUCAAIGDCEhADjIRhkGAAIEtAbseWyrftEW5V0u9YqlVBAgQIHAwgcvBxmM4BAgQIHAVkHy8uBSivNfKi5/wGgECBAjsJ2AHZD9bXyZAgMBLAo5cvcS2filKUyvru1oIECBAYJCABGQQvG4JECCwJWDXY0vlB21R3q2lXrHUKgIECBAYLHAZ3L/uCRAgQOAqIPnYaSlE+W4rO3XhswQIECDwuIAdkMetPEmAAIFdBBy52oV1/dEoTa2s72ohQIAAgU4CEpBO0LohQIDAloBdjy2VHduifLuWesVSqwgQIECgs8Clc3+6I0CAAIGrgORj0FKI0m8rg4agWwIECGQWsAOSefbFToDAEAFHroawrzuN0tTK+q4WAgQIENhJQAKyE6zPEiBAYEvArseWysC2KH3XUq9YahUBAgQI7Cxw2fn7Pk+AAAECVwHJx0GXQpRxtXLQIRoWAQIEZhKwAzLTbIqFAIFDCjhydchpWQ8qSlMr67taCBAgQOBNAhKQN0H6DAECBLYE7HpsqRy4LcrYaqlXLLWKAAECBN4scHnz93yOAAECBK4Cko+TLoUo427lpCEYNgECBI4sYAfkyLNjbAQInFLAkatTTtt60FGaWlnf1UKAAAECLwpIQF6E8xoBAgS2BOx6bKmcuC3K2GupVyy1igABAgR+KHD54fteJ0CAAIGrgORj0qUQJa5WJg1RWAQIEOgpYAekp7a+CBCYUsCRqymndR1UlKZW1ne1ECBAgMCDAhKQB6E8RoAAgS0Bux5bKhO3RYmtlnrFUqsIECBA4EmBy5PPe5wAAQIErgKSj6RLIUrcrSQlEDYBAgR+ImAH5Cd63iVAIKWAI1cpp30ddJSmVtZ3tRAgQIDAFwISkC9gNBMgQGBLwK7Hlkritiix11KvWGoVAQIECNwRuNy57zYBAgQIXAUkH5bCpkCU1lY2H9BIgAABAn8K2AH5U8NvAgQIbAg4crWBomktEKWplfVdLQQIECBwFZCAWAoECBD4RsCuxzc4bq0FojTVUq9YahUBAgQI/CVw+etvfxIgQIDAVUDyYSm8JBDlrVZe+oCXCBAgMLeAHZC551d0BAi8IODI1QtoXlkLRGlqZX1XCwECBNIKSEDSTr3ACRDYErDrsaWi7WWBKG/WUq9YahUBAgTSC1zSCwAgQIDAVUDyYSnsIhDlq63s0oGPEiBA4FwCdkDONV9GS4DADgKOXO2A6pNrgShNrazvaiFAgEAaAQlImqkWKAECWwJ2PbZUtO0mEOXLtdQrllpFgACBdAKXdBELmAABAlcByYelMEQgSq+tDBmATgkQIDBWwA7IWH+9EyAwQMCRqwHoulwLRGlqZX1XCwECBKYVkIBMO7UCI0BgS8Cux5aKtmECUXqupV6x1CoCBAhML3CZPkIBEiBA4Cog+bAUDikQZVStHHKABkWAAIH3CtgBea+nrxEgcEABR64OOCmGtBaI0tTK+q4WAgQITCMgAZlmKgVCgMCWgF2PLRVthxWIMrJa6hVLrSJAgMB0ApfpIhIQAQIErgKSD0vhlAJRRt3KKQMwaAIECHwvYAfkex93CRA4oYAjVyecNENeC0RpamV9VwsBAgROKyABOe3UGTgBAlsCdj22VLSdViDKyGupVyy1igABAqcXuJw+AgEQIEDgKiD5sBSmFIgSVStTBigoAgSyCdgByTbj4iUwoYAjVxNOqpDWAlGaWlnf1UKAAIHTCEhATjNVBkqAwJaAXY8tFW3TCkSJrJZ6xVKrCBAgcDqBy+lGbMAECBC4Ckg+LIWUAlGibiUlgKAJEDi7gB2Qs8+g8RNIKODIVcJJF/JaIEpTK+u7WggQIHBYAQnIYafGwAgQ2BKw67Gloi2tQJTIa6lXLLWKAAEChxe4HH6EBkiAAIGrgOTDUiCwIRClrZWN25oIECBwNAE7IEebEeMhQGAl4MjVikQDgbVAlKZW1ne1ECBA4DACEpDDTIWBECCwJWDXY0tFG4EvBKK011KvWGoVAQIEDidwOdyIDIgAAQJXAcmHpUDgBYEo77TywuteIUCAwN4CdkD2FvZ9AgSeFnDk6mkyLxBYC0RpamV9VwsBAgSGCUhAhtHrmACBLQG7Hlsq2gi8KBDlvVrqFUutIkCAwHCBy/ARGAABAgSuApIPS4HADgJRvtnKDp/3SQIECDwrYAfkWTHPEyDwdgFHrt5O6oME1gJRmlpZ39VCgACBbgISkG7UOiJAYEvArseWijYCOwlE+W4t9YqlVhEgQKC7wKV7jzokQIDAVUDyYSkQGCAQpc9WBnSvSwIECNgBsQYIEOgu4MhVd3IdElgLRGlqZX1XCwECBHYTkIDsRuvDBAhsCdj12FLRRmCQQJR+a6lXLLWKAAECuwtcdu9BBwQIELgKSD4sBQIHFIgyplYOODxDIkBgPgE7IPPNqYgIHE7AkavDTYkBEVgLRGlqZX1XCwECBN4mIAF5G6UPESCwJWDXY0tFG4GDCkQZVy31iqVWESBA4O0Cl7d/0QcJECBwFZB8WAoETigQZcytnHD4hkyAwPEF7IAcf46MkMDpBBy5Ot2UGTCBtUCUplbWd7UQIEDgZQEJyMt0XiRAYEvArseWijYCJxWIMu5a6hVLrSJAgMCPBS4//oIPECBA4Cog+bAUCEwoECWmViYMT0gECPQXsAPS31yPBKYTcORquikVEIG1QJSmVtZ3tRAgQOBhAQnIw1QeJEBgS8Cux5aKNgKTCkSJq5Z6xVKrCBAg8LTA5ek3vECAAIGrgOTDUiCQUCBKzK0kDF/IBAj8XMAOyM8NfYFAOgFHrtJNuYAJrAWiNLWyvquFAAECXwpIQL6kcYMAgS0Bux5bKtoIJBWIEnct9YqlVhEgQOCuwOXuEx4gQIDAVUDyYSkQILASiNLSyuqmBgIECKwF7ICsTbQQIPCXgCNXf4H4kwCBtUCUplbWd7UQIEDgt4AE5DeFHwQIbAnY9dhS0UaAwKZAlNZa6hVLrSJAgMBK4LJq0UCAAIGrgOTDUiBA4GmBKG+08vTLXiBAIIOAHZAMsyxGAk8KOHL1JJjHCRBYC0RpamV9VwsBAokFJCCJJ1/oBLYE7HpsqWgjQOAlgShv1VKvWGoVAQIEPi4MCBAg0AQkH03CvwQIvE0gypdaedtHfYgAgTML2AE58+wZO4E3CThy9SZInyFA4GuBKLda+fopdwgQSCAgAUkwyUIk8J2AXY/vdNwjQOCtAlG+Vku9YqlVBAgkFLgkjFnIBAhcBSQflgIBAt0FovTYSvfOdUiAwBEE7IAcYRaMgUBnAUeuOoPrjgCBtUCUplbWd7UQIDCxgARk4skVGoEtAbseWyraCBAYIhCl11rqFUutIkAggcAlQYxCJEDgKiD5sBQIEDicQJQRtXK4wRkQAQJ7CNgB2UPVNwkcTMCRq4NNiOEQILAWiNLUyvquFgIEJhKQgEw0mUIhsCVg12NLRRsBAocUiDKqWuoVS60iQGBCgcuEMQmJAIGrgOTDUiBA4HQCUUbcyukGb8AECDwiYAfkESXPEDiZgCNXJ5swwyVAYC0QpamV9V0tBAicWEACcuLJM3QCWwJ2PbZUtBEgcEqBKKOupV6x1CoCBCYQuEwQgxAIELgKSD4sBQIEphOIElEr0wUnIAI5BeyA5Jx3UU8m4MjVZBMqHAIE1gJRmlpZ39VCgMCJBCQgJ5osQyWwJWDXY0tFGwECUwpEiaqWesVSqwgQOKHA5YRjNmQCBK4Ckg9LgQCBdAJRIm4lXfACJjCHgB2QOeZRFMkEHLlKNuHCJUBgLRClqZX1XS0ECBxYQAJy4MkxNAJbAnY9tlS0ESCQUiBK1LXUK5ZaRYDACQQuJxijIRIgcBWQfFgKBAgQ+Esgyt+t/HXLnwQIHFPADsgx58WoCHwScOTqE4c/CBAgsBaI0tTK+q4WAgQOJCABOdBkGAqBLQG7Hlsq2ggQILAhEKWtlnrFUqsIEDigwOWAYzIkAgSuApIPS4EAAQJPCkR5vpUnX/U4AQJ9BOyA9HHWC4GnBBy5eorLwwQIEFgLRGlqZX1XCwECAwUkIAPxdU1gS8Cux5aKNgIECLwgEOWdWuoVS60iQOAAApcDjMEQCBC4Ckg+LAUCBAi8WSDK91p586d9jgCB1wTsgLzm5i0CbxVw5OqtnD5GgACBtUCUplbWd7UQINBRQALSEVtXBLYE7HpsqWgjQIDADgJRvllLvWKpVQQIDBC4DOhTlwQIXAUkH5YCAQIEOgtE6a+Vzl3rjgCBXwJ2QKwEAgMEHLkagK5LAgQI/CkQ5Y9W/mz3mwCB3QUkILsT64DAZwG7Hp89/EWAAIFhAlF6rqVesdQqAgQ6CFw69KELAgSuApIPS4EAAQIHE4gynlYONjTDITCrgB2QWWdWXIcScOTqUNNhMAQIEFgLRGlqZX1XCwECbxSQgLwR06cIbAnY9dhS0UaAAIEDCkQZUy31iqVWESCwg8Blh2/6JAECVwHJh6VAgACBkwlEGW8rJxu64RI4i4AdkLPMlHGeSsCRq1NNl8ESIEBgLRClqZX1XS0ECPxAQALyAzyvEtgSsOuxpaKNAAECJxSIMuZa6hVLrSJA4A0Clzd8wycIELgKSD4sBQIECEwmECWeViYLTTgERgnYARklr9+pBBy5mmo6BUOAAIG1QJSmVtZ3tRAg8ISABOQJLI8S2BKw67Gloo0AAQITCkSJqZZ6xVKrCBB4QeDywjteIUDgKiD5sBQIECCQTCBKvK0kC124BN4lYAfkXZK+k0rAkatU0y1YAgQIrAWiNLXy593//Prjn49/f3z8z/XGv/58wG8CBCQg1gCBJwXsejwJ5nECBAjMKhAlsFrqFaVck4/656ertktCPpH4I7eABCT3/Iv+SQHJx5NgHidAgMDsAnEN8Kvko8UvCWkS/iXwIQGxCAg8IODI1QNIHiFAgEBWgXvJR3ORhDQJ/yYXkIAkXwDCvy9g1+O+kScIECBAgAABAo8K1BOJj+btj37TcwSmEWg7H9MEJBACBAgQeLvAP/9T/gfnj17+tyCPSnluYgEJyMSTKzQCBAgQIECgg8Az/1WuBKTDhOji6AKXow/Q+AgQIECAAAECBAgQmEdAAjLPXIqEAAECBAgQGCHw6K7Go8+NiEGfBDoKSEA6YuuKAAECBAgQmFTgXnJx7/6kLMIisCUgAdlS0UaAAAECBAgQeFbgqyTjq/Znv+95ApMI1P9IPPM/nZokbGEQIECAAAECBAgQIDBCwA7ICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIAQnICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIAQnICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIAQnICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIAQnICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIAQnICHV9EiBAgAABAgQIEEgqIAFJOvHCJkCAAAECBAgQIDBCQAIyQl2fBAgQIECAAAECBJIKSECSTrywCRAgQIAAAQIECIwQkICMUNcnAQIECBAgQIAAgaQCEpCkEy9sAgQIECBAgAABAiMEJCAj1PVJgAABAgQIECBAIKmABCTpxAubAAECBAgQIECAwAgBCcgIdX0SIECAAAECBAgQSCogAUk68cImQIAAAQIECBAgMEJAAjJCXZ8ECBAgQIAAAQIEkgpIQJJOvLAJECBAgAABAgQIjBCQgIxQ1ycBAgQIECBAgACBpAISkKQTL2wCBAgQIECAAAECIwQkICPU9UmAAAECBAgQIEAgqYAEJOnEC5sAAQIECBAgQIDACAEJyAh1fRIgQIAAAQIECBBIKiABSTrxwiZAgAABAgQIECAwQkACMkJdnwQIECBAgAABAgSSCkhAkk68sAkQIECAAAECBAiMEJCAjFDXJwECBAgQIECAAIGkAhKQpBMvbAIECBAgQIAAAQIjBCQgI9T1SYAAAQIECBAgQCCpgAQk6cQLmwABAgQIECBAgMAIgf8HxRpl+Yi4HwMAAAAASUVORK5CYII=", - "text/html": [ - "" - ], - "text/plain": [ - " Size: 2MB\n", - "array([[4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " ...,\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080]], shape=(600, 800), dtype=uint32)\n", - "Coordinates:\n", - " * y (y) float64 5kB 0.003333 0.01 0.01667 0.02333 ... 3.983 3.99 3.997\n", - " * x (x) float64 6kB 0.0025 0.0075 0.0125 0.0175 ... 3.987 3.993 3.998" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ - "# find the path from start to goal,\n", - "# barriers are uncrossable cells. In this case, they are cells with a value of 0\n", - "\n", - "path_agg_8_connectivity = a_star_search(\n", - " line_agg, start, goal, barriers=[0], snap_start=True, snap_goal=True\n", + "path_agg_4 = a_star_search(\n", + " line_agg, start, goal, barriers=[0],\n", + " snap_start=True, snap_goal=True, connectivity=4\n", ")\n", "\n", - "path_shaded = dynspread(shade(path_agg_8_connectivity, cmap=[\"green\"]))\n", + "path_shaded = dynspread(shade(path_agg_4, cmap=[\"green\"]))\n", "set_background(stack(line_shaded, path_shaded, start_shaded, goal_shaded), \"black\")" ] }, @@ -243,84 +124,156 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 4-connectivity" + "### Weighted A*\n", + "\n", + "When a **friction surface** is provided, A\\* finds the *least-cost* path rather\n", + "than the geometrically shortest one. Each edge cost becomes\n", + "`geometric_distance * mean_friction_of_endpoints`, the same cost model used\n", + "by `cost_distance()`.\n", + "\n", + "This is useful for terrain routing: a path over flat ground may be longer\n", + "in distance but much cheaper than a short path over steep slopes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "For 4-connectivity distance, we use the same arguments as above, but set the connectivity to 4." + "#### Generate terrain and friction surface" ] }, { "cell_type": "code", - "execution_count": 4, - "metadata": { - "execution": { - "iopub.execute_input": "2026-01-23T01:07:50.934668Z", - "iopub.status.busy": "2026-01-23T01:07:50.934547Z", - "iopub.status.idle": "2026-01-23T01:07:51.154146Z", - "shell.execute_reply": "2026-01-23T01:07:51.153618Z" - } - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAA2ZElEQVR4Ae3dbZLjuJEG4PbGhg+8F7BPsAfeP16nPKrOZklVEotMAsiHEY6mvgjkk5gYvwOM/bdfv37969//cREgQIAAAQIECBAgQOB0gf86fQQDECBAgAABAgQIECBA4C8BAcRSIECAAAECBAgQIECgTEAAKaM2EAECBAgQIECAAAECAog1QIAAAQIECBAgQIBAmYAAUkZtIAIECBAgQIAAAQIEBBBrgAABAgQIECBAgACBMgEBpIzaQAQIECBAgAABAgQICCDWAAECBAgQIECAAAECZQICSBm1gQgQIECAAAECBAgQEECsAQIECBAgQIAAAQIEygQEkDJqAxEgQIAAAQIECBAgIIBYAwQIECBAgAABAgQIlAkIIGXUBiJAgAABAgQIECBAQACxBggQIECAAAECBAgQKBMQQMqoDUSAAAECBAgQIECAgABiDRAgQIAAAQIECBAgUCYggJRRG4gAAQIECBAgQIAAAQHEGiBAgAABAgQIECBAoExAACmjNhABAgQIECBAgAABAgKINUCAAAECBAgQIECAQJmAAFJGbSACBAgQIECAAAECBAQQa4AAAQIECBAgQIAAgTIBAaSM2kAECBAgQIAAAQIECAgg1gABAgQIECBAgAABAmUCAkgZtYEIECBAgAABAgQIEBBArAECBAgQIECAAAECBMoEBJAyagMRIECAAAECBAgQICCAWAMECBAgQIAAAQIECJQJCCBl1AYiQIAAAQIECBAgQEAAsQYIECBAgAABAgQIECgTEEDKqA1EgAABAgQIECBAgIAAYg0QIECAAAECBAgQIFAmIICUURuIAAECBAgQIECAAAEBxBogQIAAAQIECBAgQKBMQAApozYQAQIECBAgQIAAAQICiDVAgAABAgQIECBAgECZgABSRm0gAgQIECBAgAABAgQEEGuAAAECBAgQIECAAIEyAQGkjNpABAgQIECAAAECBAgIINYAAQIECBAgQIAAAQJlAgJIGbWBCBAgQIAAAQIECBAQQKwBAgQIECBAgAABAgTKBASQMmoDESBAgAABAgQIECAggFgDBAgQIECAAAECBAiUCQggZdQGIkCAAAECBAgQIEBAALEGCBAgQIAAAQIECBAoExBAyqgNRIAAAQIECBAgQICAAGINECBAgAABAgQIECBQJiCAlFEbiAABAgQIECBAgAABAcQaIECAAAECBAgQIECgTEAAKaM2EAECBAgQIECAAAECAog1QIAAAQIECBAgQIBAmYAAUkZtIAIECBAgQIAAAQIEBBBrgAABAgQIECBAgACBMgEBpIzaQAQIECBAgAABAgQICCDWAAECBAgQIECAAAECZQICSBm1gQgQIECAAAECBAgQEECsAQIECBAgQIAAAQIEygQEkDJqAxEgQIAAAQIECBAgIIBYAwQIECBAgAABAgQIlAkIIGXUBiJAgAABAgQIECBAQACxBggQIECAAAECBAgQKBMQQMqoDUSAAAECBAgQIECAgABiDRAgQIAAAQIECBAgUCYggJRRG4gAAQIECBAgQIAAAQHEGiBAgAABAgQIECBAoExAACmjNhABAgQIECBAgAABAgKINUCAAAECBAgQIECAQJmAAFJGbSACBAgQIECAAAECBAQQa4AAAQIECBAgQIAAgTIBAaSM2kAECBAgQIAAAQIECAgg1gABAgQIECBAgAABAmUCAkgZtYEIECBAgAABAgQIEBBArAECBAgQIECAAAECBMoEBJAyagMRIECAAAECBAgQICCAWAMECBAgQIAAAQIECJQJCCBl1AYiQIAAAQIECBAgQEAAsQYIECBAgAABAgQIECgTEEDKqA1EgAABAgQIECBAgIAAYg0QIECAAAECBAgQIFAmIICUURuIAAECBAgQIECAAAEBxBogQIAAAQIECBAgQKBMQAApozYQAQIECBAgQIAAAQICiDVAgAABAgQIECBAgECZgABSRm0gAgQIECBAgAABAgQEEGuAAAECBAgQIECAAIEyAQGkjNpABAgQIECAAAECBAgIINYAAQIECBAgQIAAAQJlAgJIGbWBCBAgQIAAAQIECBAQQKwBAgQIECBAgAABAgTKBASQMmoDESBAgAABAgQIECAggFgDBAgQIECAAAECBAiUCQggZdQGIkCAAAECBAgQIEBAALEGCBAgQIAAAQIECBAoExBAyqgNRIAAAQIECBAgQICAAGINECBAgAABAgQIECBQJiCAlFEbiAABAgQIECBAgAABAcQaIECAAAECBAgQIECgTEAAKaM2EAECBAgQIECAAAECAog1QIAAAQIECBAgQIBAmYAAUkZtIAIECBAgQIAAAQIEBBBrgAABAgQIECBAgACBMgEBpIzaQAQIECBAgAABAgQICCDWAAECBAgQIECAAAECZQICSBm1gQgQIECAAAECBAgQEECsAQIECBAgQIAAAQIEygQEkDJqAxEgQIAAAQIECBAgIIBYAwQIECBAgAABAgQIlAkIIGXUBiJAgAABAgQIECBAQACxBggQIECAAAECBAgQKBMQQMqoDUSAAAECBAgQIECAgABiDRAgQIAAAQIECBAgUCYggJRRG4gAAQIECBAgQIAAAQHEGiBAgAABAgQIECBAoExAACmjNhABAgQIECBAgAABAgKINUCAAAECBAgQIECAQJmAAFJGbSACBAgQIECAAAECBAQQa4AAAQIECBAgQIAAgTIBAaSM2kAECBAgQIAAAQIECAgg1gABAgQIECBAgAABAmUCAkgZtYEIECBAgAABAgQIEBBArAECBAgQIECAAAECBMoEBJAyagMRIECAAAECBAgQICCAWAMECBAgQIAAAQIECJQJCCBl1AYiQIAAAQIECBAgQEAAsQYIECBAgAABAgQIECgTEEDKqA1EgAABAgQIECBAgIAAYg0QIECAAAECBAgQIFAmIICUURuIAAECBAgQIECAAAEBxBogQIAAAQIECBAgQKBMQAApozYQAQIECBAgQIAAAQICiDVAgAABAgQIECBAgECZgABSRm0gAgQIECBAgAABAgQEEGuAAAECBAgQIECAAIEyAQGkjNpABAgQIECAAAECBAgIINYAAQIECBAgQIAAAQJlAgJIGbWBCBAgQIAAAQIECBAQQKwBAgQIECBAgAABAgTKBASQMmoDESBAgAABAgQIECAggFgDBAgQIECAAAECBAiUCQggZdQGIkCAAAECBAgQIEBAALEGCBAgQIAAAQIECBAoExBAyqgNRIAAAQIECBAgQIDAfyMgQOCzwL8+v/Xrbw/e8xYBAgQIECBAgMB7AnZA3vPy7QYCj8JHlP3s/QYkSiRAgAABAgQIHCYggBxG6UErCHwXMr77fAUDNRAgQIAAAQIEzhQQQM7U9eypBF4NF69+b6riTZYAAQIECBAgUCQggBRBG4YAAQIECBAgQIAAgV+//EvoVgGBHQL/98//2fErPyFAgACBTgJ//+f/dipXrQReFhBAXqbyRQK/BfxN5beFu74CEcT9tdC3/yp/LOAfUD128S6BLOAIVtZwT4AAAQIECBDYKXAP5YL5TkA/ayMggLRptUK/E3j1/+fj1e99N57PCRAgQGAdgXv4WKcilRA4T0AAOc/WkycU+C5cxOe21ydsrCkTIEDgJIH4e4LwcRKuxy4r4N8BWba1CtsrECHj0f/U7j2cxNb6PYTYZt+r7HcECBCYX0DwmL+HKrhGwA7INe5GHVwgwsbtP//868/NfCN45CCy+dhLAgQIEFhcQPhYvMHKO1VAADmV18NXFxBCVu+w+ggQIPCnQAQP4eNPE68IvCvgCNa7Yr5PYCOQQ0jcuwgQIEBgTQHBY82+qqpewA5IvbkRFxSI4JGDyIIlKokAAQKtBYSP1u1X/MECAsjBoB7XW0AI6d1/1RMgsJ6AI1fr9VRF1ws4gnV9D8xgMYEcQuLeRYAAAQJzCtj1mLNvZj2+gB2Q8XtkhhMKRPDIQWTCEkyZAAECrQWEj9btV/zJAgLIycAe31tACOndf9UTIDCfgCNX8/XMjOcTcARrvp6Z8WQCOYTEvYsAAQIExhSw6zFmX8xqPQE7IOv1VEUDCkTwyEFkwCmaEgECBFoLCB+t26/4YgEBpBjccL0FhJDe/Vc9AQLjCThyNV5PzGh9AUew1u+xCgcTyCEk7l0ECBAgcI2AXY9r3I1KwA6INUDgAoEIHjmIXDAFQxIgQKC1gPDRuv2Kv1hAALm4AYbvLSCE9O6/6gkQqBdw5Kre3IgEtgKOYG1FvCZQLJBDSNy7CBAgQOAcAbse57h6KoF3BeyAvCvm+wROEIjgkYPICUN4JAECBFoLCB+t26/4wQQEkMEaYjq9BYSQ3v1XPQECxws4cnW8qScS+KmAI1g/FfR7AgcL5BAS9y4CBAgQ2Cdg12Ofm18ROFvADsjZwp5PYIdABI8cRHY8wk8IECDQWkD4aN1+xQ8uIIAM3iDT6y0ghPTuv+oJEHhfwJGr9838gkC1gCNY1eLGI/CmQA4hce8iQIAAgccCdj0eu3iXwGgCdkBG64j5EHggEMEjB5EHX/EWAQIEWgsIH63br/jJBASQyRpmur0FhJDe/Vc9AQKfBRy5+mziHQKjCziCNXqHzI/ARiCHkLh3ESBAoKuAXY+unVf37AJ2QGbvoPm3FIjgkYNISwRFEyDQWkD4aN1+xU8uIIBM3kDT7y0ghPTuv+oJdBRw5Kpj19W8moAjWKt1VD3tBHIIiXsXAQIEVhWw67FqZ9XVTcAOSLeOq3dJgQgeOYgsWaSiCBBoLSB8tG6/4hcTEEAWa6hyegsIIb37r3oCKwo4crViV9XUXcARrO4rQP3LCeQQEvcuAgQIzCpg12PWzpk3ga8F7IB87eNTAlMKRPDIQWTKIkyaAIHWAsJH6/YrfnEBAWTxBiuvt4AQ0rv/qicwo4AjVzN2zZwJvCfgCNZ7Xr5NYDqBHELi3kWAAIFRBex6jNoZ8yJwrIAdkGM9PY3AkAIRPHIQGXKSJkWAQGsB4aN1+xXfTEAAadZw5fYWEEJ691/1BEYUcORqxK6YE4FzBRzBOtfX0wkMJ5BDSNy7CBAgcJWAXY+r5I1L4FoBOyDX+hudwCUCETxyELlkEgYlQKC1gPDRuv2Kby4ggDRfAMrvLSCE9O6/6glcIeDI1RXqxiQwloAjWGP1w2wIlAvkEBL3LgIECJwlYNfjLFnPJTCXgB2QufpltgROEYjgkYPIKYN4KAECrQWEj9btVzyBPwQEkD84vCDQW0AI6d1/1RM4Q8CRqzNUPZPA3AKOYM3dP7MncLhADiFx7yJAgMBeAbsee+X8jsDaAnZA1u6v6gjsEojgkYPIrof4EQECrQWEj9btVzyBLwUEkC95fEigt4AQ0rv/qiewR8CRqz1qfkOgl4AjWL36rVoCbwvkEBL3LgIECDwTsOvxTMb7BAhkATsgWcM9AQIPBSJ45CDy8EveJECgtYDw0br9iifwloAA8haXLxPoLSCE9O6/6gk8EnDk6pGK9wgQ+ErAEayvdHxGgMAngRxC4t5FgEBfAbsefXuvcgI/EbAD8hM9vyXQVCCCRw4iTRmUTaC1gPDRuv2KJ/AjAQHkR3x+TKC3gBDSu/+q7yngyFXPvquawJECjmAdqelZBBoK5BAS9y4CBNYVsOuxbm9VRqBSwA5IpbaxCCwqEMEjB5FFy1QWgdYCwkfr9iuewKECAsihnB5GoLeAENK7/6pfU8CRqzX7qioCVwo4gnWlvrEJLCiQQ0jcuwgQmFfArse8vTNzAiML2AEZuTvmRmBSgQgeOYhMWoZpE2gtIHy0br/iCZwqIICcyuvhBHoLCCG9+6/6OQUcuZqzb2ZNYCYBR7Bm6pa5EphQIIeQuHcRIDCugF2PcXtjZgRWErADslI31UJgUIEIHjmIDDpN0yLQWkD4aN1+xRMoFRBASrkNRqC3gBDSu/+qH1PAkasx+2JWBFYWcARr5e6qjcCAAjmExL2LAIHrBOx6XGdvZAKdBeyAdO6+2glcJBDBIweRi6ZhWAKtBYSP1u1XPIFLBQSQS/kNTqC3gBDSu/+qv0bAkatr3I1KgMBvAUewflu4I0DgAoEcQuLeRYDAeQJ2Pc6z9WQCBF4XsAPyupVvEiBwkkAEjxxEThrGYwm0FhA+Wrdf8QSGEhBAhmqHyRDoLSCE9O6/6s8RcOTqHFdPJUBgv4AjWPvt/JIAgRMEcgiJexcBAvsF7Hrst/NLAgTOE7ADcp6tJxMgsFMggkcOIjsf42cEWgsIH63br3gCQwsIIEO3x+QI9BYQQnr3X/X7BBy52ufmVwQI1Ak4glVnbSQCBHYI5BAS9y4CBJ4L2PV4buMTAgTGEbADMk4vzIQAgScCETxyEHnyNW8TaC0gfLRuv+IJTCUggEzVLpMl0FtACOndf9U/FnDk6rGLdwkQGFfAEaxxe2NmBAg8EMghJO5dBDoL2PXo3H21E5hXwA7IvL0zcwJtBSJ45CDSFkLhrQWEj9btVzyBqQUEkKnbZ/IEegsIIb3737V6R666dl7dBNYRcARrnV6qhEBLgRxC4t5FYGUBux4rd1dtBPoI2AHp02uVElhWIIJHDiLLFqqw1gLCR+v2K57AUgICyFLtVAyB3gJCSO/+r1q9I1erdlZdBPoKOILVt/cqJ7CkQA4hce8iMLOAXY+Zu2fuBAg8E7AD8kzG+wQITCsQwSMHkWkLMfHWAsJH6/YrnsDSAgLI0u1VHIHeAkJI7/7PWr0jV7N2zrwJEHhVwBGsV6V8jwCBKQVyCIl7F4GRBex6jNwdcyNA4CgBOyBHSXoOAQLDCkTwyEFk2ImaWGsB4aN1+xVPoJWAANKq3Yol0FtACOnd/1Grd+Rq1M6YFwECZwk4gnWWrOcSIDCkQA4hce8icKWAXY8r9Y1NgMBVAnZArpI3LgEClwlE8MhB5LKJGLi1gPDRuv2KJ9BaQABp3X7FE+gtIIT07v9V1TtydZW8cQkQGEXAEaxROmEeBAhcIpBDSNy7CJwpYNfjTF3PJkBgFgE7ILN0yjwJEDhNIIJHDiKnDeTBrQWEj9btVzwBAklAAEkYbgkQ6C0ghPTu/1nVO3J1lqznEiAwq4AjWLN2zrwJEDhFIIeQuHcR+ImAXY+f6PktAQKrCtgBWbWz6iJAYLdABI8cRHY/yA9bCwgfrduveAIEvhAQQL7A8REBAr0FhJDe/d9bvSNXe+X8jgCBLgKOYHXptDoJENglkENI3LsIfCVg1+MrHZ8RIEDgPwJ2QKwEAgQIfCMQwSMHkW++7uOmAsJH08YrmwCBtwUEkLfJ/IAAga4CQkjXzn9dtyNXX/v4lAABAlsBR7C2Il4TIEDgC4EcQuLe1VvArkfv/queAIF9AnZA9rn5FQECjQUieOQg0piidenCR+v2K54AgR8ICCA/wPNTAgR6CwghPfvvyFXPvquaAIHjBBzBOs7SkwgQaCiQQ0jcu9YWsOuxdn9VR4BAjYAdkBpnoxAgsLBABI8cRBYutXVpwkfr9iueAIEDBQSQAzE9igCB3gJCyJr9d+Rqzb6qigCB6wQcwbrO3sgECCwokENI3LvmFrDrMXf/zJ4AgTEF7ICM2RezIkBgYoEIHjmITFxK66kLH63br3gCBE4UEEBOxPVoAgR6Cwghc/bfkas5+2bWBAjMI+AI1jy9MlMCBCYUyCEk7l1jC9j1GLs/ZkeAwBoCdkDW6KMqCBAYWCCCRw4iA0+19dSEj9btVzwBAoUCAkghtqEIEOgtIISM2X9Hrsbsi1kRILCugCNY6/ZWZQQIDCiQQ0jcu64VsOtxrb/RCRDoKWAHpGffVU2AwIUCETxyELlwKq2HFj5at1/xBAhcKCCAXIhvaAIEegsIIdf035Gra9yNSoAAgbuAI1h3CX8SIEDgAoEcQuLeda6AXY9zfT2dAAECrwjYAXlFyXcIECBwokAEjxxEThyq9aOFj9btVzwBAgMJCCADNcNUCBDoLSCEnNN/R67OcfVUAgQI7BVwBGuvnN8RIEDgBIEcQuLe9TMBux4/8/NrAgQInCFgB+QMVc8kQIDADwQieOQg8oNHtf6p8NG6/YonQGBgAQFk4OaYGgECvQWEkH39d+Rqn5tfESBAoErAEawqaeMQIEBgh0AOIXHv+lrArsfXPj4lQIDACAJ2QEbogjkQIEDgC4EIHjmIfPHV1h8JH63br3gCBCYSEEAmapapEiDQW0AIedx/R64eu3iXAAECowo4gjVqZ8yLAAECDwRyCIn77pddj+4rQP0ECMwoYAdkxq6ZMwECrQUieOQg0hVD+OjaeXUTIDC7gAAyewfNnwCBtgJdQ4gjV22XvMIJEFhEwBGsRRqpDAIEegrkEBL3q192PVbvsPoIEOggYAekQ5fVSIDA0gIRPHIQWbVY4WPVzqqLAIFuAgJIt46rlwCBZQVWDSGOXC27ZBVGgEBTAUewmjZe2QQIrCmQQ0jcz37Z9Zi9g+ZPgACBzwJ2QD6beIcAAQJTC0TwyEFk1mKEj1k7Z94ECBD4WkAA+drHpwQIEJhWYNYQ4sjVtEvOxAkQIPCSgCNYLzH5EgECBOYUyCEk7ke/7HqM3iHzI0CAwM8F7ID83NATCBAgMLRABI8cREadrPAxamfMiwABAscKCCDHenoaAQIEhhUYNYQ4cjXskjExAgQInCLgCNYprB5KgACBMQVyCIn7qy+7Hld3wPgECBCoF7ADUm9uRAIECFwqEMEjB5GrJiN8XCVvXAIECFwrIIBc6290AgQIXCZwVQhx5OqylhuYAAECQwg4gjVEG0yCAAEC1wjkEBL3Z192Pc4W9nwCBAiML2AHZPwemSEBAgROFYjgkYPIWYMJH2fJei4BAgTmEhBA5uqX2RIgQOA0gbNCiCNXp7XMgwkQIDClgCNYU7bNpAkQIHCOQA4hcf/Ty67HTwX9ngABAusJ2AFZr6cqIkCAwI8EInjkILL3YcLHXjm/I0CAwNoCAsja/VUdAQIEdgvsDSGOXO0m90MCBAi0EHAEq0WbFUmAAIF9AjmExP13l12P74R8ToAAAQJ2QKwBAgQIEPhSIIJHDiLPvix8PJPxPgECBAhkAQEka7gnQIAAgacCz0KII1dPyXxAgAABAg8EHMF6gOItAgQIEHgssA0hdj0eO3mXAAECBJ4LCCDPbXxC4CYQ/wXLRYDAcwF/jTy38QkBAgQIfBYQQD6beIfAHwLxT3xdBAj8FtgGDn+N/LZxR+AusP3r5P6+PwkQ+PXLvwNiFRAgQIDAywLxX6oicNxDR/wZ7/kvWy8T+iIBAgTaCwgg7ZcAAAIECLwmcA8f22/fA4kQspXxmgABAgQeCQggj1S8R4AAAQIfAvcdjvuux8cHm5v7bsjmbS8JECBAgMAfAv4dkD84vCBAgACBLPBs1yN/J9/nEPJdYMm/c0+AAAECfQTsgPTptUoJECDwlsC74eP+8AgeOYjc3/cnAQIECBAIAQHEOiBAgACBPwRePXL1x48evBBCHqB4iwABAgR+OYJlERAgQIDAh8DeXY+PB2xucgiJexcBAgQIELADYg0QIECAwE3g6PBxZ43gkYPI/X1/EiBAgEBPAQGkZ99VTYAAgQ+Bo45cfTzwyY0Q8gTG2wQIEGgm4AhWs4YrlwABAlngrF2PPEa+zyEk7l0ECBAg0E/ADki/nquYAAECN4Hq8HFnj+CRg8j9fX8SIECAQA8BAaRHn1VJgACBD4GqI1cfAz65EUKewHibAAECiws4grV4g5VHgACBLHDVrkeeQ77PISTuXQQIECCwvoAdkPV7rEICBAjcBEYLH/e2RPDIQeT+vj8JECBAYE0BAWTNvqqKAAECHwKjHLn6mNCTGyHkCYy3CRAgsJiAI1iLNVQ5BAgQyAKj7nrkOeb7HELi3kWAAAEC6wnYAVmvpyoiQIDATWC28HFvWwSPHETu7/uTAAECBNYQEEDW6KMqCBAg8CEwy5Grjwk/uRFCnsB4mwABApMLOII1eQNNnwABAllg1l2PXEO+zyEk7l0ECBAgML+AHZD5e6gCAgQI3ARWCx/3tkbwyEHk/r4/CRAgQGBOAQFkzr6ZNQECBD4EVjly9VHQkxsh5AmMtwkQIDCZgCNYkzXMdAkQIJAFVt31yDXm+xxC4t5FgAABAvMJ2AGZr2dmTIAAgZtAt/Bxb3sEjxxE7u/7kwABAgTmEBBA5uiTWRIgQOBDoMuRq4+Cn9wIIU9gvE2AAIHBBRzBGrxBpkeAAIEs0HXXIxvk+xxC4t5FgAABAuML2AEZv0dmSIAAgZuA8PF4IUTwyEHk8be8S4AAAQKjCAggo3TCPAgQIPBEwJGrJzCbt4WQDYiXBAgQGFTAEaxBG2NaBAgQCAG7Hu+tgxxC4t5FgAABAuMJ2AEZrydmRIAAgZuA8LFvIUTwyEFk31P8igABAgTOEhBAzpL1XAIECOwUcORqJ9zmZ0LIBsRLAgQIDCLgCNYgjTANAgQIhIBdj2PXQQ4hce8iQIAAgesF7IBc3wMzIECAwE1A+DhnIUTwyEHknFE8lQABAgReFRBAXpXyPQIECJwk4MjVSbCbxwohGxAvCRAgcJGAI1gXwRuWAAECIWDXo3Yd5BAS9y4CBAgQqBewA1JvbkQCBAjcBISPaxZCBI8cRK6ZhVEJECDQV0AA6dt7lRMgcJGAI1cXwW+GFUI2IF4SIECgSMARrCJowxAgQCAE7HqMtQ5yCIl7FwECBAicL2AH5HxjIxAgQOAmIHyMuRAieOQgMuYszYoAAQLrCAgg6/RSJQQIDCrgyNWgjdlMSwjZgHhJgACBkwQcwToJ1mMJECAQAnY95loHOYTEvYsAAQIEjhewA3K8qScSIEDgJiB8zLkQInjkIDJnFWZNgACBcQUEkHF7Y2YECEwq4MjVpI3bTFsI2YB4SYAAgYMEHME6CNJjCBAgEAJ2PdZaBzmExL2LAAECBH4uYAfk54aeQIAAgZuA8LHmQojgkYPImlWqigABAnUCAkidtZEIEFhUwJGrRRu7KUsI2YB4SYAAgZ0CjmDthPMzAgQIhIBdj17rIIeQuHcRIECAwPsCdkDeN/MLAgQI3ASEj54LIYJHDiI9FVRNgACB/QICyH47vyRAoKmAI1dNG78pWwjZgHhJgACBFwUcwXoRytcIECAQAnY9rIMskENI3LsIECBA4HsBOyDfG/kGAQIEbgLCh4XwSCCCRw4ij77jPQIECBD4LSCA/LZwR4AAgYcCjlw9ZPHmRkAI2YB4SYAAgScCjmA9gfE2AQIEQsCuh3XwjkAOIXHvIkCAAIHPAnZAPpt4hwABAjcB4cNC2CMQwSMHkT3P8BsCBAisLCCArNxdtREgsEvAkatdbH60ERBCNiBeEiBA4C8BR7AsBQIECCQBux4Jw+2PBXIIiXsXAQIECPz6ZQfEKiBAgMBfAsKHpXCGQASPHETOGMMzCRAgMJOAADJTt8yVAIFTBBy5OoXVQzcCQsgGxEsCBNoKOILVtvUKJ0AgBOx6WAeVAjmExL2LAAECHQXsgHTsupoJELgJCB8WwhUCETxyELliDsYkQIDAlQICyJX6xiZA4BIBR64uYTfoRkAI2YB4SYBAGwFHsNq0WqEECISAXQ/rYCSBHELi3kWAAIEOAnZAOnRZjQQI3ASEDwthRIEIHjmIjDhHcyJAgMCRAgLIkZqeRYDAkAKOXA3ZFpPaCAghGxAvCRBYVsARrGVbqzACBELArod1MJNADiFx7yJAgMCKAnZAVuyqmggQuAkIHxbCjAIRPHIQmbEGcyZAgMBXAgLIVzo+I0BgSgFHrqZsm0lvBISQDYiXBAgsI+AI1jKtVAgBAiFg18M6WEkgh5C4dxEgQGAFATsgK3RRDQQI3ASEDwthRYEIHjmIrFijmggQ6CUggPTqt2oJLCngyNWSbVXURkAI2YB4SYDAtAKOYE3bOhMnQCAE7HpYB50EcgiJexcBAgRmFLADMmPXzJkAgZuA8GEhdBSI4JGDSEcDNRMgMLeAADJ3/8yeQEsBR65atl3RGwEhZAPiJQEC0wg4gjVNq0yUAIEQsOthHRD4LZBDSNy7CBAgMIOAHZAZumSOBAjcBIQPC4HAZ4EIHjmIfP6GdwgQIDCWgAAyVj/MhgCBBwKOXD1A8RaBjYAQsgHxkgCBYQUcwRq2NSZGgEAI2PWwDgi8LpBDSNy7CBAgMKKAHZARu2JOBAjcBIQPC4HA+wIRPHIQef8JfkGAAIFzBQSQc309nQCBHQKOXO1A8xMCGwEhZAPiJQECwwg4gjVMK0yEAIEQsOthHRA4TiCHkLh3ESBAYAQBOyAjdMEcCBC4CQgfFgKB4wUieOQgcvwInkiAAIH3BASQ97x8mwCBEwQcuToB1SMJbASEkA2IlwQIXCbgCNZl9AYmQCAE7HpYBwTqBHIIiXsXAQIErhCwA3KFujEJELgJCB8WAoF6gQgeOYjUz8CIBAh0FxBAuq8A9RO4QMCRqwvQDUlgIyCEbEC8JECgTMARrDJqAxEgEAJ2PawDAuMI5BAS9y4CBAhUCNgBqVA2BgECNwHhw0IgMJ5ABI8cRMaboRkRILCagACyWkfVQ2BAAUeuBmyKKRHYCAghGxAvCRA4TcARrNNoPZgAgRCw62EdEJhHIIeQuHcRIEDgDAE7IGeoeiYBAjcB4cNCIDCfQASPHETmq8CMCRAYXUAAGb1D5kdgQgFHriZsmikT2AgIIRsQLwkQOEzAEazDKD2IAIEQsOthHRBYRyCHkLh3ESBA4AgBOyBHKHoGAQI3AeHDQiCwnkAEjxxE1qtQRQQIVAsIINXixiOwoIAjVws2VUkENgJCyAbESwIEdgs4grWbzg8JEAgBux7WAYE+AjmExL2LAAECewTsgOxR8xsCBG4CwoeFQKCfQASPHET6CaiYAIGfCgggPxX0ewINBRy5ath0JRPYCAghGxAvCRB4WcARrJepfJEAgRCw62EdECBwF8ghJO5dBAgQeEXADsgrSr5DgMBNQPiwEAgQ2ApE8MhBZPu51wQIENgKCCBbEa8JEPgk4MjVJxJvECCwERBCNiBeEiDwVMARrKc0PiBAIATselgHBAi8KpBDSNy7CBAg8EjADsgjFe8RIHATED4sBAIE3hWI4JGDyLu/930CBNYXEEDW77EKCbwt4MjV22R+QIDARkAI2YB4SYDAh4AjWB8UbggQCAG7HtYBAQJHCeQQEvcuAgQIhIAdEOuAAIEPAeHjg8INAQIHCUTwyEHkoMd6DAECEwsIIBM3z9QJHCXgyNVRkp5DgMAzASHkmYz3CfQTcASrX89VTOAPAbsef3B4QYDAiQI5hMS9iwCBngJ2QHr2XdUEbgLCh4VAgEC1QASPHESqxzceAQLXCwgg1/fADAiUCzhyVU5uQAIENgJCyAbESwKNBBzBatRspRIIAbse1gEBAqMI5BAS9y4CBHoI2AHp0WdVErgJCB8WAgECowlE8MhBZLT5mQ8BAscLCCDHm3oigeEEHLkariUmRIDARkAI2YB4SWBhAUewFm6u0giEgF0P64AAgVkEcgiJexcBAmsK2AFZs6+qInATED4sBAIEZhOI4JGDyGzzN18CBL4XEEC+N/INAtMJOHI1XctMmACBjYAQsgHxksBCAo5gLdRMpRAIAbse1gEBAqsI5BAS9y4CBNYQsAOyRh9VQeAmIHxYCAQIrCYQwSMHkdXqUw+BjgICSMeuq3k5AUeulmupgggQ2AgIIRsQLwlMLOAI1sTNM3UCIWDXwzogQKCLQA4hce8iQGBOATsgc/bNrAncBIQPC4EAgW4CETxyEOlWv3oJrCAggKzQRTW0E3Dkql3LFUyAwEZACNmAeElgIgFHsCZqlqkSCAG7HtYBAQIE/iOQQ0jcuwgQmEPADsgcfTJLAjcB4cNCIECAwJ8CETxyEPnzU68IEBhRQAAZsSvmRGAj4MjVBsRLAgQIbASEkA2IlwQGFnAEa+DmmBqBELDrYR0QIEDgNYEcQuLeRYDAmAJ2QMbsi1kRuAkIHxYCAQIE3hOI4JGDyHu/9m0CBCoEBJAKZWMQeFPAkas3wXydAAECGwEhZAPiJYGBBBzBGqgZpkIgBOx6WAcECBA4RiCHkLh3ESAwhoAdkDH6YBYEbgLCh4VAgACBYwUieOQgcuzTPY0AgT0CAsgeNb8hcLCAI1cHg3ocAQIENgJCyAbESwIXCjiCdSG+oQmEgF0P64AAAQI1AjmExL2LAIFrBOyAXONuVAI3AeHDQiBAgECtQASPHERqRzcaAQIhIIBYBwQuEHDk6gJ0QxIgQCAJCCEJwy2BYgFHsIrBDUfAroc1QIAAgTEEcgiJexcBAjUCdkBqnI1C4CYgfFgIBAgQGEsggkcOImPNzmwIrCkggKzZV1UNJuDI1WANMR0CBAhsBISQDYiXBE4UcATrRFyPJhACdj2sAwIECMwhkENI3LsIEDhHwA7IOa6eSuAmIHxYCAQIEJhLIIJHDiJzzd5sCcwhIIDM0SeznEzAkavJGma6BAgQ2AgIIRsQLwkcKOAI1oGYHkUgBOx6WAcECBBYQyCHkLh3ESBwjIAdkGMcPYXATUD4sBAIECCwlkAEjxxE1qpONQSuERBArnE36mICjlwt1lDlECBAYCMghGxAvCTwAwFHsH6A56cEQsCuh3VAgACBHgI5hMS9iwCBfQJ2QPa5+RWBm4DwYSEQIECgl0AEjxxEelWvWgLHCNgBOcbRU5oJRPCIyz8Ba9Z45RIgQOAvgXsI+fT3gX/95wt///XvHZJ//PXlv2EjQCALCCBZwz2BFwTseryA5CsECBBoIHAPIVHqLYj8FT4+lR7vCyGfWLzRV0AA6dt7le8QED52oPkJAQIEFhb42AF5Fj7utQshdwl/EvglgFgEBF4QcOTqBSRfIUCAQFeB78LH3UUIuUv4s7mAANJ8ASj/ewG7Ht8b+QYBAgQIECBA4FWBOJH4am5/9Zm+R2AdgX+uU4pKCBAgQOAkgX+88Vz/LsgbWL66qoAAsmpn1UWAAAECBAjUCLzzj3IFkJqeGGVoAf8/IEO3x+QIECBAgAABAgQIrCUggKzVT9UQIECAAAEC1QKv7mq8+r3q+RuPQLGAAFIMbjgCBAgQIEBgQYHvwsV3ny9IoiQCzwQEkGcy3idAgAABAgQIvCPwLGQ8e/+dZ/sugYUE4i+Jd/7VqYVKVwoBAgQIECBAgAABAtUCdkCqxY1HgAABAgQIECBAoLGAANK4+UonQIAAAQIECBAgUC0ggFSLG48AAQIECBAgQIBAYwEBpHHzlU6AAAECBAgQIECgWkAAqRY3HgECBAgQIECAAIHGAgJI4+YrnQABAgQIECBAgEC1gABSLW48AgQIECBAgAABAo0FBJDGzVc6AQIECBAgQIAAgWoBAaRa3HgECBAgQIAAAQIEGgsIII2br3QCBAgQIECAAAEC1QICSLW48QgQIECAAAECBAg0FhBAGjdf6QQIECBAgAABAgSqBQSQanHjESBAgAABAgQIEGgsIIA0br7SCRAgQIAAAQIECFQLCCDV4sYjQIAAAQIECBAg0FhAAGncfKUTIECAAAECBAgQqBYQQKrFjUeAAAECBAgQIECgsYAA0rj5SidAgAABAgQIECBQLSCAVIsbjwABAgQIECBAgEBjAQGkcfOVToAAAQIECBAgQKBaQACpFjceAQIECBAgQIAAgcYCAkjj5iudAAECBAgQIECAQLWAAFItbjwCBAgQIECAAAECjQUEkMbNVzoBAgQIECBAgACBagEBpFrceAQIECBAgAABAgQaCwggjZuvdAIECBAgQIAAAQLVAgJItbjxCBAgQIAAAQIECDQWEEAaN1/pBAgQIECAAAECBKoFBJBqceMRIECAAAECBAgQaCwggDRuvtIJECBAgAABAgQIVAsIINXixiNAgAABAgQIECDQWEAAadx8pRMgQIAAAQIECBCoFhBAqsWNR4AAAQIECBAgQKCxgADSuPlKJ0CAAAECBAgQIFAtIIBUixuPAAECBAgQIECAQGMBAaRx85VOgAABAgQIECBAoFpAAKkWNx4BAgQIECBAgACBxgICSOPmK50AAQIECBAgQIBAtYAAUi1uPAIECBAgQIAAAQKNBQSQxs1XOgECBAgQIECAAIFqAQGkWtx4BAgQIECAAAECBBoLCCCNm690AgQIECBAgAABAtUCAki1uPEIECBAgAABAgQINBYQQBo3X+kECBAgQIAAAQIEqgUEkGpx4xEgQIAAAQIECBBoLCCANG6+0gkQIECAAAECBAhUCwgg1eLGI0CAAAECBAgQINBYQABp3HylEyBAgAABAgQIEKgWEECqxY1HgAABAgQIECBAoLGAANK4+UonQIAAAQIECBAgUC0ggFSLG48AAQIECBAgQIBAYwEBpHHzlU6AAAECBAgQIECgWkAAqRY3HgECBAgQIECAAIHGAgJI4+YrnQABAgQIECBAgEC1gABSLW48AgQIECBAgAABAo0FBJDGzVc6AQIECBAgQIAAgWoBAaRa3HgECBAgQIAAAQIEGgsIII2br3QCBAgQIECAAAEC1QICSLW48QgQIECAAAECBAg0FhBAGjdf6QQIECBAgAABAgSqBQSQanHjESBAgAABAgQIEGgsIIA0br7SCRAgQIAAAQIECFQLCCDV4sYjQIAAAQIECBAg0FhAAGncfKUTIECAAAECBAgQqBYQQKrFjUeAAAECBAgQIECgsYAA0rj5SidAgAABAgQIECBQLSCAVIsbjwABAgQIECBAgEBjAQGkcfOVToAAAQIECBAgQKBaQACpFjceAQIECBAgQIAAgcYCAkjj5iudAAECBAgQIECAQLWAAFItbjwCBAgQIECAAAECjQUEkMbNVzoBAgQIECBAgACBagEBpFrceAQIECBAgAABAgQaCwggjZuvdAIECBAgQIAAAQLVAgJItbjxCBAgQIAAAQIECDQWEEAaN1/pBAgQIECAAAECBKoFBJBqceMRIECAAAECBAgQaCwggDRuvtIJECBAgAABAgQIVAsIINXixiNAgAABAgQIECDQWEAAadx8pRMgQIAAAQIECBCoFhBAqsWNR4AAAQIECBAgQKCxgADSuPlKJ0CAAAECBAgQIFAtIIBUixuPAAECBAgQIECAQGMBAaRx85VOgAABAgQIECBAoFpAAKkWNx4BAgQIECBAgACBxgICSOPmK50AAQIECBAgQIBAtYAAUi1uPAIECBAgQIAAAQKNBQSQxs1XOgECBAgQIECAAIFqAQGkWtx4BAgQIECAAAECBBoLCCCNm690AgQIECBAgAABAtUCAki1uPEIECBAgAABAgQINBYQQBo3X+kECBAgQIAAAQIEqgUEkGpx4xEgQIAAAQIECBBoLCCANG6+0gkQIECAAAECBAhUCwgg1eLGI0CAAAECBAgQINBYQABp3HylEyBAgAABAgQIEKgWEECqxY1HgAABAgQIECBAoLGAANK4+UonQIAAAQIECBAgUC0ggFSLG48AAQIECBAgQIBAYwEBpHHzlU6AAAECBAgQIECgWkAAqRY3HgECBAgQIECAAIHGAgJI4+YrnQABAgQIECBAgEC1gABSLW48AgQIECBAgAABAo0FBJDGzVc6AQIECBAgQIAAgWoBAaRa3HgECBAgQIAAAQIEGgsIII2br3QCBAgQIECAAAEC1QICSLW48QgQIECAAAECBAg0FhBAGjdf6QQIECBAgAABAgSqBQSQanHjESBAgAABAgQIEGgsIIA0br7SCRAgQIAAAQIECFQLCCDV4sYjQIAAAQIECBAg0FhAAGncfKUTIECAAAECBAgQqBYQQKrFjUeAAAECBAgQIECgsYAA0rj5SidAgAABAgQIECBQLSCAVIsbjwABAgQIECBAgEBjAQGkcfOVToAAAQIECBAgQKBaQACpFjceAQIECBAgQIAAgcYCAkjj5iudAAECBAgQIECAQLWAAFItbjwCBAgQIECAAAECjQUEkMbNVzoBAgQIECBAgACBagEBpFrceAQIECBAgAABAgQaCwggjZuvdAIECBAgQIAAAQLVAgJItbjxCBAgQIAAAQIECDQWEEAaN1/pBAgQIECAAAECBKoFBJBqceMRIECAAAECBAgQaCwggDRuvtIJECBAgAABAgQIVAsIINXixiNAgAABAgQIECDQWEAAadx8pRMgQIAAAQIECBCoFhBAqsWNR4AAAQIECBAgQKCxgADSuPlKJ0CAAAECBAgQIFAtIIBUixuPAAECBAgQIECAQGMBAaRx85VOgAABAgQIECBAoFpAAKkWNx4BAgQIECBAgACBxgICSOPmK50AAQIECBAgQIBAtYAAUi1uPAIECBAgQIAAAQKNBQSQxs1XOgECBAgQIECAAIFqAQGkWtx4BAgQIECAAAECBBoLCCCNm690AgQIECBAgAABAtUCAki1uPEIECBAgAABAgQINBYQQBo3X+kECBAgQIAAAQIEqgUEkGpx4xEgQIAAAQIECBBoLCCANG6+0gkQIECAAAECBAhUCwgg1eLGI0CAAAECBAgQINBYQABp3HylEyBAgAABAgQIEKgW+H85950PbG3GKQAAAABJRU5ErkJggg==", - "text/html": [ - "" - ], - "text/plain": [ - " Size: 2MB\n", - "array([[4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " ...,\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080],\n", - " [4278190080, 4278190080, 4278190080, ..., 4278190080, 4278190080,\n", - " 4278190080]], shape=(600, 800), dtype=uint32)\n", - "Coordinates:\n", - " * y (y) float64 5kB 0.003333 0.01 0.01667 0.02333 ... 3.983 3.99 3.997\n", - " * x (x) float64 6kB 0.0025 0.0075 0.0125 0.0175 ... 3.987 3.993 3.998" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "# Generate a terrain raster\nW, H = 600, 400\ntemplate = xr.DataArray(np.zeros((H, W)), dims=['y', 'x'])\nterrain = generate_terrain(template, seed=12345)\n\nshade(terrain, cmap=Elevation, how=\"linear\")" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "# Derive friction from slope: steeper = higher cost\nterrain_slope = slope(terrain)\n\n# slope() produces NaN at edges; fill with 0 so friction = 1 there (passable, flat)\n# Add 1 so flat areas have friction=1 rather than 0\nfriction = terrain_slope.fillna(0.0) + 1.0\n\nshade(friction, cmap=[\"lightgreen\", \"yellow\", \"red\"], how=\"eq_hist\")" + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "# find the path from start to goal,\n", - "# barriers are uncrossable cells. In this case, they are cells with a value of 0\n", + "#### Find paths: unweighted vs weighted\n", "\n", - "path_agg_4_connectivity = a_star_search(\n", - " line_agg, start, goal, barriers=[0], snap_start=True, snap_goal=True, connectivity=4\n", + "We pick a start and goal on opposite sides of a ridge, then compare\n", + "the path found with and without friction weighting." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose start and goal coordinates\n", + "start_coord = (400.0, 50.0)\n", + "goal_coord = (100.0, 450.0)\n", + "\n", + "# Unweighted path (geometric distance only)\n", + "path_unweighted = a_star_search(terrain, start_coord, goal_coord)\n", + "\n", + "# Weighted path (friction from slope)\n", + "path_weighted = a_star_search(terrain, start_coord, goal_coord, friction=friction)\n", + "\n", + "print(f\"Unweighted path cost at goal: {np.nanmax(path_unweighted):.2f}\")\n", + "print(f\"Weighted path cost at goal: {np.nanmax(path_weighted):.2f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Visualise both paths on the terrain\n", + "\n", + "Red = unweighted (geometrically shortest), green = weighted (least-cost).\n", + "Notice how the weighted path avoids steep ridges even if it's longer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "terrain_shaded = shade(terrain, cmap=Elevation, how=\"linear\", alpha=180)\n", + "\n", + "unweighted_shaded = dynspread(\n", + " shade(path_unweighted, cmap=[\"red\", \"red\"], how=\"linear\", min_alpha=255),\n", + " threshold=1, max_px=2,\n", + ")\n", + "weighted_shaded = dynspread(\n", + " shade(path_weighted, cmap=[\"lime\", \"lime\"], how=\"linear\", min_alpha=255),\n", + " threshold=1, max_px=2,\n", ")\n", "\n", - "path_shaded = dynspread(shade(path_agg_4_connectivity, cmap=[\"green\"]))\n", - "set_background(stack(line_shaded, path_shaded, start_shaded, goal_shaded), \"black\")" + "stack(terrain_shaded, unweighted_shaded, weighted_shaded)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### References" + "#### Validate against cost_distance\n", + "\n", + "For a single source, the A\\* path cost at the goal should match the\n", + "`cost_distance()` value at that same pixel. This confirms both algorithms\n", + "use the same cost model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a source raster with a single target at the start pixel\n", + "source_raster = terrain.copy(data=np.zeros_like(terrain.data))\n", + "# Find the start pixel index\n", + "from xrspatial.pathfinding import _get_pixel_id\n", + "spy, spx = _get_pixel_id(start_coord, terrain)\n", + "source_raster.data[spy, spx] = 1.0\n", + "\n", + "cd = cost_distance(source_raster, friction)\n", + "\n", + "# Compare: A* cost at goal vs cost_distance at goal\n", + "gpy, gpx = _get_pixel_id(goal_coord, terrain)\n", + "cd_val = float(cd.values[gpy, gpx])\n", + "astar_val = float(path_weighted.values[gpy, gpx])\n", + "\n", + "print(f\"cost_distance at goal: {cd_val:.4f}\")\n", + "print(f\"A* path cost at goal: {astar_val:.4f}\")\n", + "print(f\"Match: {np.isclose(cd_val, astar_val, rtol=1e-5)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "- A* search algorithm: https://en.wikipedia.org/wiki/A*_search_algorithm" + "### References\n", + "\n", + "- A\\* search algorithm: https://en.wikipedia.org/wiki/A*_search_algorithm\n", + "- Red Blob Games β€” A\\* implementation: https://www.redblobgames.com/pathfinding/a-star/implementation.html\n", + "- Cost distance (weighted proximity): `xrspatial.cost_distance`" ] + }, + { + "cell_type": "markdown", + "source": "### Dask (out-of-core)\n\nFor rasters that don't fit in memory, wrap your data in a Dask array.\nA\\* explores only the pixels it needs (typically O(path_length) to\nO(path_length^2)), loading chunks on demand through an LRU cache.\nThe output is a lazy Dask array of the same chunk structure.", + "metadata": {} + }, + { + "cell_type": "code", + "source": "import dask.array as da\n\n# Re-use the terrain and friction from above, but as dask arrays\nterrain_dask = terrain.copy()\nterrain_dask.data = da.from_array(terrain.data, chunks=(200, 300))\n\nfriction_dask = friction.copy()\nfriction_dask.data = da.from_array(friction.data, chunks=(200, 300))\n\n# Run A* β€” chunks are loaded on demand, not all at once\npath_dask = a_star_search(terrain_dask, start_coord, goal_coord, friction=friction_dask)\nprint(f\"Result type: {type(path_dask.data)}\")\nprint(f\"Chunks: {path_dask.data.chunks}\")\n\n# Compute to verify it matches the in-memory result\ndask_cost = float(path_dask.values[gpy, gpx])\nprint(f\"Dask A* cost at goal: {dask_cost:.4f}\")\nprint(f\"NumPy A* cost at goal: {astar_val:.4f}\")\nprint(f\"Match: {np.isclose(dask_cost, astar_val, rtol=1e-5)}\")", + "metadata": {}, + "execution_count": null, + "outputs": [] } ], "metadata": { @@ -344,4 +297,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/xrspatial/pathfinding.py b/xrspatial/pathfinding.py index 8c271bca..7d393734 100644 --- a/xrspatial/pathfinding.py +++ b/xrspatial/pathfinding.py @@ -1,10 +1,24 @@ +import heapq import warnings +from collections import OrderedDict +from math import sqrt from typing import Optional, Union import numpy as np import xarray as xr -from xrspatial.utils import get_dataarray_resolution, ngjit +try: + import dask.array as da + import dask +except ImportError: + da = None + dask = None + +from xrspatial.cost_distance import _heap_push, _heap_pop +from xrspatial.utils import ( + get_dataarray_resolution, ngjit, + has_cuda_and_cupy, is_cupy_array, is_dask_cupy, has_dask_array, +) NONE = -1 @@ -42,6 +56,16 @@ def _is_not_crossable(cell_value, barriers): return False +def _is_not_crossable_py(cell_value, barriers): + """Pure Python version of _is_not_crossable for the dask A* loop.""" + if np.isnan(cell_value): + return True + for b in barriers: + if cell_value == b: + return True + return False + + @ngjit def _is_inside(py, px, h, w): inside = True @@ -52,36 +76,6 @@ def _is_inside(py, px, h, w): return inside -@ngjit -def _distance(x1, y1, x2, y2): - # euclidean distance in pixel space from (y1, x1) to (y2, x2) - return np.sqrt((x1 - x2)**2 + (y1 - y2)**2) - - -@ngjit -def _heuristic(x1, y1, x2, y2): - # heuristic to estimate distance between 2 point - # TODO: what if we want to use another distance metric? - return _distance(x1, y1, x2, y2) - - -@ngjit -def _min_cost_pixel_id(cost, is_open): - height, width = cost.shape - py = NONE - px = NONE - # set min cost to a very big number - # this value is only an estimation - min_cost = (height + width) ** 2 - for i in range(height): - for j in range(width): - if is_open[i, j] and cost[i, j] < min_cost: - min_cost = cost[i, j] - py = i - px = j - return py, px - - @ngjit def _find_nearest_pixel(py, px, data, barriers): # if the cell is already valid, return itself @@ -89,15 +83,15 @@ def _find_nearest_pixel(py, px, data, barriers): return py, px height, width = data.shape - # init min distance as max possible distance - min_distance = _distance(0, 0, height - 1, width - 1) + # init min distance as max possible distance (pixel space) + min_distance = np.sqrt(float((height - 1) ** 2 + (width - 1) ** 2)) # return of the function nearest_y = NONE nearest_x = NONE for y in range(height): for x in range(width): if not _is_not_crossable(data[y, x], barriers): - d = _distance(x, y, px, py) + d = np.sqrt(float((x - px) ** 2 + (y - py) ** 2)) if d < min_distance: min_distance = d nearest_y = y @@ -106,11 +100,34 @@ def _find_nearest_pixel(py, px, data, barriers): return nearest_y, nearest_x +def _neighborhood_structure(cellsize_x, cellsize_y, connectivity=8): + """Return (dy, dx, dd) with cellsize-scaled geometric distances.""" + if connectivity == 8: + dy = np.array([-1, -1, -1, 0, 0, 1, 1, 1], dtype=np.int64) + dx = np.array([-1, 0, 1, -1, 1, -1, 0, 1], dtype=np.int64) + dd = np.array([ + sqrt(cellsize_y ** 2 + cellsize_x ** 2), # (-1,-1) + cellsize_y, # (-1, 0) + sqrt(cellsize_y ** 2 + cellsize_x ** 2), # (-1,+1) + cellsize_x, # ( 0,-1) + cellsize_x, # ( 0,+1) + sqrt(cellsize_y ** 2 + cellsize_x ** 2), # (+1,-1) + cellsize_y, # (+1, 0) + sqrt(cellsize_y ** 2 + cellsize_x ** 2), # (+1,+1) + ], dtype=np.float64) + else: + dy = np.array([0, -1, 1, 0], dtype=np.int64) + dx = np.array([-1, 0, 0, 1], dtype=np.int64) + dd = np.array([cellsize_x, cellsize_y, cellsize_y, cellsize_x], + dtype=np.float64) + return dy, dx, dd + + @ngjit -def _reconstruct_path(path_img, parent_ys, parent_xs, cost, +def _reconstruct_path(path_img, parent_ys, parent_xs, g_cost, start_py, start_px, goal_py, goal_px): # construct path output image as a 2d array with NaNs for non-path pixels, - # and the value of the path pixels being the current cost up to that point + # and the value of the path pixels being the g-cost up to that point current_x = goal_px current_y = goal_py @@ -118,11 +135,11 @@ def _reconstruct_path(path_img, parent_ys, parent_xs, cost, parent_ys[current_y, current_x] != NONE: # exist path from start to goal # add cost at start - path_img[start_py, start_px] = cost[start_py, start_px] + path_img[start_py, start_px] = g_cost[start_py, start_px] # add cost along the path while current_x != start_px or current_y != start_py: # value of a path pixel is the cost up to that point - path_img[current_y, current_x] = cost[current_y, current_x] + path_img[current_y, current_x] = g_cost[current_y, current_x] parent_y = parent_ys[current_y, current_x] parent_x = parent_xs[current_y, current_x] current_y = parent_y @@ -130,26 +147,16 @@ def _reconstruct_path(path_img, parent_ys, parent_xs, cost, return -def _neighborhood_structure(connectivity=8): - if connectivity == 8: - # 8-connectivity - neighbor_xs = [-1, -1, -1, 0, 0, 1, 1, 1] - neighbor_ys = [-1, 0, 1, -1, 1, -1, 0, 1] - else: - # 4-connectivity - neighbor_ys = [0, -1, 1, 0] - neighbor_xs = [-1, 0, 0, 1] - return np.array(neighbor_ys), np.array(neighbor_xs) - - @ngjit def _a_star_search(data, path_img, start_py, start_px, goal_py, goal_px, - barriers, neighbor_ys, neighbor_xs): + barriers, dy, dx, dd, friction, f_min, use_friction, + cellsize_x, cellsize_y): height, width = data.shape + n_neighbors = len(dy) + # parent of the (i, j) pixel is the pixel at # (parent_ys[i, j], parent_xs[i, j]) - # first initialize parent of all cells as invalid (NONE, NONE) parent_ys = np.ones((height, width), dtype=np.int64) * NONE parent_xs = np.ones((height, width), dtype=np.int64) * NONE @@ -157,77 +164,328 @@ def _a_star_search(data, path_img, start_py, start_px, goal_py, goal_px, parent_ys[start_py, start_px] = start_py parent_xs[start_py, start_px] = start_px - # distance from start to the current node - d_from_start = np.zeros_like(data, dtype=np.float64) - # total cost of the node: cost = d_from_start + d_to_goal - # heuristic β€” estimated distance from the current node to the end node - cost = np.zeros_like(data, dtype=np.float64) + # g-cost: distance from start to the current node + g_cost = np.full((height, width), np.inf, dtype=np.float64) - # initialize both open and closed list all False - is_open = np.zeros(data.shape, dtype=np.bool_) - is_closed = np.zeros(data.shape, dtype=np.bool_) + visited = np.zeros((height, width), dtype=np.int8) + + # Heap arrays + max_heap = height * width + h_keys = np.empty(max_heap, dtype=np.float64) + h_rows = np.empty(max_heap, dtype=np.int64) + h_cols = np.empty(max_heap, dtype=np.int64) + h_size = 0 if not _is_not_crossable(data[start_py, start_px], barriers): - # if start node is crossable - # add the start node to open list - is_open[start_py, start_px] = True - # init cost at start location - d_from_start[start_py, start_px] = 0 - cost[start_py, start_px] = d_from_start[start_py, start_px] + \ - _heuristic(start_px, start_py, goal_px, goal_py) - - num_open = np.sum(is_open) - while num_open > 0: - py, px = _min_cost_pixel_id(cost, is_open) - # pop current node off open list, add it to closed list - is_open[py][px] = 0 - is_closed[py][px] = True + # Check friction at start when using friction + if use_friction: + f_start_val = friction[start_py, start_px] + if not (np.isfinite(f_start_val) and f_start_val > 0.0): + return + + g_cost[start_py, start_px] = 0.0 + + # Compute heuristic for start + dy_goal = abs(start_py - goal_py) * cellsize_y + dx_goal = abs(start_px - goal_px) * cellsize_x + h = np.sqrt(dy_goal ** 2 + dx_goal ** 2) + if use_friction: + h *= f_min + + h_size = _heap_push(h_keys, h_rows, h_cols, h_size, + h, start_py, start_px) + + while h_size > 0: + f_u, py, px, h_size = _heap_pop(h_keys, h_rows, h_cols, h_size) + + if visited[py, px]: + continue + visited[py, px] = 1 + # found the goal - if (py, px) == (goal_py, goal_px): - # reconstruct path + if py == goal_py and px == goal_px: _reconstruct_path(path_img, parent_ys, parent_xs, - d_from_start, start_py, start_px, + g_cost, start_py, start_px, goal_py, goal_px) return + g_u = g_cost[py, px] + # visit neighborhood - for y, x in zip(neighbor_ys, neighbor_xs): - neighbor_y = py + y - neighbor_x = px + x + for i in range(n_neighbors): + ny = py + dy[i] + nx = px + dx[i] - # neighbor is within the surface image - if neighbor_y > height - 1 or neighbor_y < 0 \ - or neighbor_x > width - 1 or neighbor_x < 0: + if ny < 0 or ny >= height or nx < 0 or nx >= width: continue - - # walkable - if _is_not_crossable(data[neighbor_y][neighbor_x], barriers): + if visited[ny, nx]: continue + if _is_not_crossable(data[ny, nx], barriers): + continue + + # Compute edge cost + if use_friction: + f_u_val = friction[py, px] + f_v_val = friction[ny, nx] + # impassable if friction is NaN or non-positive + if not (np.isfinite(f_v_val) and f_v_val > 0.0): + continue + edge_cost = dd[i] * (f_u_val + f_v_val) * 0.5 + else: + edge_cost = dd[i] + + new_g = g_u + edge_cost + + if new_g < g_cost[ny, nx]: + g_cost[ny, nx] = new_g + parent_ys[ny, nx] = py + parent_xs[ny, nx] = px + + # Compute heuristic + dy_goal = abs(ny - goal_py) * cellsize_y + dx_goal = abs(nx - goal_px) * cellsize_x + h = np.sqrt(dy_goal ** 2 + dx_goal ** 2) + if use_friction: + h *= f_min + + f_val = new_g + h + h_size = _heap_push(h_keys, h_rows, h_cols, h_size, + f_val, ny, nx) - # check if neighbor is in the closed list - if is_closed[neighbor_y, neighbor_x]: + return + + +# --------------------------------------------------------------------------- +# LRU chunk cache for dask A* +# --------------------------------------------------------------------------- + +class _ChunkCache: + """OrderedDict-based LRU cache for dask chunks.""" + + def __init__(self, maxsize=128): + self._cache = OrderedDict() + self._maxsize = maxsize + + def get(self, key, loader): + """Return cached chunk or call *loader()*, evicting oldest if full.""" + if key in self._cache: + self._cache.move_to_end(key) + return self._cache[key] + value = loader() + if len(self._cache) >= self._maxsize: + self._cache.popitem(last=False) + self._cache[key] = value + return value + + +# --------------------------------------------------------------------------- +# Sparse dask A* +# --------------------------------------------------------------------------- + +def _a_star_dask(surface_da, friction_da, start_py, start_px, + goal_py, goal_px, barriers, dy, dx, dd, + f_min, use_friction, cellsize_x, cellsize_y, is_cupy): + """Run A* on a dask-backed array, loading chunks on demand. + + Returns a list of (row, col, cost) tuples for path pixels, + or [] if no path exists. + """ + height, width = surface_da.shape + n_neighbors = len(dy) + + # Chunk boundaries (cumulative sums of chunk sizes) + row_chunks = np.array(surface_da.chunks[0]) + col_chunks = np.array(surface_da.chunks[1]) + row_bounds = np.cumsum(row_chunks) + col_bounds = np.cumsum(col_chunks) + + surface_cache = _ChunkCache() + friction_cache = _ChunkCache() if use_friction else None + + def _load_chunk(da_arr, cache, iy, ix): + """Load and cache a single chunk, converting cupy->numpy.""" + def loader(): + block = da_arr.blocks[iy, ix].compute() + if is_cupy: + block = block.get() + return np.asarray(block, dtype=np.float64) + return cache.get((iy, ix), loader) + + def _get_value(da_arr, cache, r, c): + """Get a scalar value at global (r, c) via chunk cache.""" + iy = int(np.searchsorted(row_bounds, r, side='right')) + ix = int(np.searchsorted(col_bounds, c, side='right')) + chunk = _load_chunk(da_arr, cache, iy, ix) + local_r = r - (int(row_bounds[iy - 1]) if iy > 0 else 0) + local_c = c - (int(col_bounds[ix - 1]) if ix > 0 else 0) + return float(chunk[local_r, local_c]) + + # Check start + start_val = _get_value(surface_da, surface_cache, start_py, start_px) + if _is_not_crossable_py(start_val, barriers): + return [] + + if use_friction: + f_start = _get_value(friction_da, friction_cache, start_py, start_px) + if not (np.isfinite(f_start) and f_start > 0.0): + return [] + + # A* data structures (sparse β€” dict/set, not full arrays) + g_cost = {(start_py, start_px): 0.0} + parent = {} + visited = set() + + counter = 0 # tie-breaker for stable heap ordering + + # Heuristic for start + dy_goal = abs(start_py - goal_py) * cellsize_y + dx_goal = abs(start_px - goal_px) * cellsize_x + h = sqrt(dy_goal ** 2 + dx_goal ** 2) + if use_friction: + h *= f_min + + heap = [(h, counter, start_py, start_px)] + + while heap: + f_u, _, py, px = heapq.heappop(heap) + + if (py, px) in visited: + continue + visited.add((py, px)) + + # Found goal β€” reconstruct path + if py == goal_py and px == goal_px: + path = [] + cr, cc = goal_py, goal_px + path.append((cr, cc, g_cost[(cr, cc)])) + while (cr, cc) in parent: + cr, cc = parent[(cr, cc)] + path.append((cr, cc, g_cost[(cr, cc)])) + return path + + g_u = g_cost[(py, px)] + + for i in range(n_neighbors): + ny = py + int(dy[i]) + nx = px + int(dx[i]) + + if ny < 0 or ny >= height or nx < 0 or nx >= width: + continue + if (ny, nx) in visited: continue - # distance from start to this neighbor - d = d_from_start[py, px] + _distance(px, py, - neighbor_x, neighbor_y) - # if neighbor is already in the open list - if is_open[neighbor_y, neighbor_x] and \ - d > d_from_start[neighbor_y, neighbor_x]: + n_val = _get_value(surface_da, surface_cache, ny, nx) + if _is_not_crossable_py(n_val, barriers): continue - # calculate cost - d_from_start[neighbor_y, neighbor_x] = d - d_to_goal = _heuristic(neighbor_x, neighbor_y, goal_px, goal_py) - cost[neighbor_y, neighbor_x] = \ - d_from_start[neighbor_y, neighbor_x] + d_to_goal - # add neighbor to the open list - is_open[neighbor_y, neighbor_x] = True - parent_ys[neighbor_y, neighbor_x] = py - parent_xs[neighbor_y, neighbor_x] = px - - num_open = np.sum(is_open) - return + if use_friction: + f_u_val = _get_value(friction_da, friction_cache, py, px) + f_v_val = _get_value(friction_da, friction_cache, ny, nx) + if not (np.isfinite(f_v_val) and f_v_val > 0.0): + continue + edge_cost = float(dd[i]) * (f_u_val + f_v_val) * 0.5 + else: + edge_cost = float(dd[i]) + + new_g = g_u + edge_cost + + if new_g < g_cost.get((ny, nx), float('inf')): + g_cost[(ny, nx)] = new_g + parent[(ny, nx)] = (py, px) + + dy_goal = abs(ny - goal_py) * cellsize_y + dx_goal = abs(nx - goal_px) * cellsize_x + h = sqrt(dy_goal ** 2 + dx_goal ** 2) + if use_friction: + h *= f_min + + counter += 1 + heapq.heappush(heap, (new_g + h, counter, ny, nx)) + + return [] # no path + + +# --------------------------------------------------------------------------- +# Sparse path β†’ lazy dask output +# --------------------------------------------------------------------------- + +def _path_to_dask_array(path_pixels, shape, chunks, is_cupy): + """Convert sparse path list to a lazy dask array of the original shape. + + *path_pixels* is a list of ``(row, col, cost)`` tuples. + Non-path pixels are NaN. + """ + height, width = shape + row_chunks = chunks[0] + col_chunks = chunks[1] + row_bounds = np.cumsum(row_chunks) + col_bounds = np.cumsum(col_chunks) + + # Group path pixels by chunk + chunk_pixels = {} # {(iy, ix): [(local_r, local_c, cost), ...]} + for r, c, cost in path_pixels: + iy = int(np.searchsorted(row_bounds, r, side='right')) + ix = int(np.searchsorted(col_bounds, c, side='right')) + local_r = r - (int(row_bounds[iy - 1]) if iy > 0 else 0) + local_c = c - (int(col_bounds[ix - 1]) if ix > 0 else 0) + chunk_pixels.setdefault((iy, ix), []).append( + (local_r, local_c, cost)) + + n_row_chunks = len(row_chunks) + n_col_chunks = len(col_chunks) + + if is_cupy: + import cupy + + @dask.delayed + def _make_block_cupy(ch, cw, pixels): + block = np.full((ch, cw), np.nan, dtype=np.float64) + for lr, lc, cost in pixels: + block[lr, lc] = cost + return cupy.asarray(block) + + blocks = [] + for iy in range(n_row_chunks): + row = [] + for ix in range(n_col_chunks): + ch = int(row_chunks[iy]) + cw = int(col_chunks[ix]) + pixels = chunk_pixels.get((iy, ix), []) + row.append( + da.from_delayed( + _make_block_cupy(ch, cw, pixels), + shape=(ch, cw), + dtype=np.float64, + meta=cupy.array((), dtype=np.float64), + ) + ) + blocks.append(row) + else: + @dask.delayed + def _make_block(ch, cw, pixels): + block = np.full((ch, cw), np.nan, dtype=np.float64) + for lr, lc, cost in pixels: + block[lr, lc] = cost + return block + + blocks = [] + for iy in range(n_row_chunks): + row = [] + for ix in range(n_col_chunks): + ch = int(row_chunks[iy]) + cw = int(col_chunks[ix]) + pixels = chunk_pixels.get((iy, ix), []) + row.append( + da.from_delayed( + _make_block(ch, cw, pixels), + shape=(ch, cw), + dtype=np.float64, + meta=np.array((), dtype=np.float64), + ) + ) + blocks.append(row) + + return da.block(blocks) def a_star_search(surface: xr.DataArray, @@ -238,21 +496,39 @@ def a_star_search(surface: xr.DataArray, y: Optional[str] = 'y', connectivity: int = 8, snap_start: bool = False, - snap_goal: bool = False) -> xr.DataArray: + snap_goal: bool = False, + friction: xr.DataArray = None) -> xr.DataArray: """ - Calculate distance from a starting point to a goal through a - surface graph. Starting location and goal location should be within - the graph. + Calculate the least-cost path from a starting point to a goal through + a surface graph, optionally weighted by a friction surface. + + A* is a modification of Dijkstra's Algorithm that is optimized for + a single destination. It prioritizes paths that seem to be leading + closer to a goal using an admissible heuristic. + + When a friction surface is provided, edge costs are + ``geometric_distance * mean_friction_of_endpoints``, matching the + cost model used by :func:`cost_distance`. The heuristic is scaled + by the minimum friction value to remain admissible. + + The output is an equal-sized ``xr.DataArray`` with NaN for non-path + pixels and the accumulated cost at each path pixel. - A* is a modification of Dijkstra’s Algorithm that is optimized for - a single destination. Dijkstra’s Algorithm can find paths to all - locations; A* finds paths to one location, or the closest of several - locations. It prioritizes paths that seem to be leading closer to - a goal. + **Backend support** - The output is an equal sized Xarray.DataArray with NaNs for non-path - pixels, and the value of the path pixels being the current cost up - to that point. + ============= =========================================================== + Backend Strategy + ============= =========================================================== + NumPy Numba-jitted kernel (fast, in-memory) + Dask Sparse Python A* with LRU chunk cache β€” loads chunks on + demand so the full grid never needs to fit in RAM + CuPy CPU fallback (transfers to numpy, runs numba kernel, + transfers back) + Dask + CuPy Same sparse A* as Dask, with cupyβ†’numpy chunk conversion + ============= =========================================================== + + ``snap_start`` and ``snap_goal`` are not supported with Dask-backed + arrays (raises ``ValueError``). Parameters ---------- @@ -267,7 +543,7 @@ def a_star_search(surface: xr.DataArray, (cannot cross). x : str, default='x' Name of the x coordinate in input surface raster. - y: str, default='x' + y: str, default='y' Name of the y coordinate in input surface raster. connectivity : int, default=8 snap_start: bool, default=False @@ -276,6 +552,12 @@ def a_star_search(surface: xr.DataArray, snap_goal: bool, default=False Snap the goal location to the nearest valid value before beginning pathfinding. + friction : xr.DataArray, optional + 2-D friction (cost) surface. Must have the same shape as + *surface*. Values must be positive and finite for passable + cells; NaN or ``<= 0`` marks impassable barriers. When + provided, edge costs become + ``geometric_distance * mean_friction_of_endpoints``. Returns ------- @@ -312,16 +594,6 @@ def a_star_search(surface: xr.DataArray, >>> start = (3, 0) >>> goal = (0, 1) >>> path_agg = a_star_search(agg, start, goal, barriers, 'lon', 'lat') - >>> print(path_agg) - - array([[ nan, nan, nan, nan], - [0. , nan, nan, nan], - [ nan, 1.41421356, nan, nan], - [ nan, nan, 2.82842712, nan], - [ nan, 4.24264069, nan, nan]]) - Coordinates: - * lon (lon) float64 0.0 1.0 2.0 3.0 - * lat (lat) float64 4.0 3.0 2.0 1.0 0.0 """ if surface.ndim != 2: @@ -334,6 +606,21 @@ def a_star_search(surface: xr.DataArray, if connectivity != 4 and connectivity != 8: raise ValueError("Use either 4 or 8-connectivity.") + # Detect backend + surface_data = surface.data + _is_dask = da is not None and isinstance(surface_data, da.Array) + _is_cupy_backend = ( + not _is_dask + and has_cuda_and_cupy() + and is_cupy_array(surface_data) + ) + _is_dask_cupy = _is_dask and has_cuda_and_cupy() and is_dask_cupy(surface) + + # compute cellsize + cellsize_x, cellsize_y = get_dataarray_resolution(surface, x, y) + cellsize_x = abs(float(cellsize_x)) + cellsize_y = abs(float(cellsize_y)) + # convert starting and ending point from geo coords to pixel coords start_py, start_px = _get_pixel_id(start, surface, x, y) goal_py, goal_px = _get_pixel_id(goal, surface, x, y) @@ -348,33 +635,142 @@ def a_star_search(surface: xr.DataArray, barriers = np.array(barriers) - if snap_start: - # find nearest valid pixel to the start location - start_py, start_px = _find_nearest_pixel( - start_py, start_px, surface.data, barriers - ) - if _is_not_crossable(surface.data[start_py, start_px], barriers): - warnings.warn("Start at a non crossable location", Warning) - - if snap_goal: - # find nearest valid pixel to the goal location - goal_py, goal_px = _find_nearest_pixel( - goal_py, goal_px, surface.data, barriers + # --- Snap / crossability checks --- + if _is_dask: + # Snapping requires O(h*w) scan β€” not supported for dask + if snap_start: + raise ValueError( + "snap_start is not supported with dask-backed arrays; " + "ensure the start pixel is valid before calling a_star_search" + ) + if snap_goal: + raise ValueError( + "snap_goal is not supported with dask-backed arrays; " + "ensure the goal pixel is valid before calling a_star_search" + ) + # Single-pixel crossability check via .compute() + start_val = float(surface_data[start_py, start_px].compute()) + if _is_not_crossable_py(start_val, barriers): + warnings.warn("Start at a non crossable location", Warning) + goal_val = float(surface_data[goal_py, goal_px].compute()) + if _is_not_crossable_py(goal_val, barriers): + warnings.warn("End at a non crossable location", Warning) + elif _is_cupy_backend: + # CuPy: use .get() for scalar access in numpy-land + surface_np = surface_data.get() + if snap_start: + start_py, start_px = _find_nearest_pixel( + start_py, start_px, surface_np, barriers + ) + if _is_not_crossable(surface_np[start_py, start_px], barriers): + warnings.warn("Start at a non crossable location", Warning) + if snap_goal: + goal_py, goal_px = _find_nearest_pixel( + goal_py, goal_px, surface_np, barriers + ) + if _is_not_crossable(surface_np[goal_py, goal_px], barriers): + warnings.warn("End at a non crossable location", Warning) + else: + # numpy path + if snap_start: + start_py, start_px = _find_nearest_pixel( + start_py, start_px, surface_data, barriers + ) + if _is_not_crossable(surface_data[start_py, start_px], barriers): + warnings.warn("Start at a non crossable location", Warning) + if snap_goal: + goal_py, goal_px = _find_nearest_pixel( + goal_py, goal_px, surface_data, barriers + ) + if _is_not_crossable(surface_data[goal_py, goal_px], barriers): + warnings.warn("End at a non crossable location", Warning) + + # Build neighborhood with cellsize-scaled distances + dy, dx, dd = _neighborhood_structure(cellsize_x, cellsize_y, connectivity) + + # --- Handle friction --- + if friction is not None: + if friction.shape != surface.shape: + raise ValueError("friction must have the same shape as surface") + use_friction = True + else: + use_friction = False + + # --- Backend dispatch --- + if _is_dask: + # Dask or dask+cupy path + if use_friction: + friction_data = friction.data + # Rechunk friction to match surface if needed + if isinstance(friction_data, da.Array): + friction_data = friction_data.rechunk(surface_data.chunks) + else: + friction_data = da.from_array( + friction_data, chunks=surface_data.chunks) + # f_min via dask (same pattern as cost_distance) + positive_friction = da.where( + friction_data > 0, friction_data, np.inf) + f_min = float(da.nanmin(positive_friction).compute()) + if not (np.isfinite(f_min) and f_min > 0): + raise ValueError("friction has no positive finite values") + else: + friction_data = None + f_min = 1.0 + + path_pixels = _a_star_dask( + surface_data, friction_data, + start_py, start_px, goal_py, goal_px, + barriers, dy, dx, dd, + f_min, use_friction, cellsize_x, cellsize_y, + _is_dask_cupy, ) - if _is_not_crossable(surface.data[goal_py, goal_px], barriers): - warnings.warn("End at a non crossable location", Warning) - - # 2d output image that stores the path - path_img = np.zeros_like(surface, dtype=np.float64) - # first, initialize all cells as np.nans - path_img[:] = np.nan + path_data = _path_to_dask_array( + path_pixels, surface.shape, surface_data.chunks, _is_dask_cupy) + + elif _is_cupy_backend: + import cupy + # Transfer to CPU, run numpy kernel, transfer back + if 'surface_np' not in dir(): + surface_np = surface_data.get() + if use_friction: + friction_np = np.asarray(friction.data.get(), dtype=np.float64) + mask = np.isfinite(friction_np) & (friction_np > 0) + if not np.any(mask): + raise ValueError("friction has no positive finite values") + f_min = float(np.min(friction_np[mask])) + else: + friction_np = np.ones((h, w), dtype=np.float64) + f_min = 1.0 + + path_img = np.full(surface.shape, np.nan, dtype=np.float64) + if start_py != NONE: + _a_star_search(surface_np, path_img, start_py, start_px, + goal_py, goal_px, barriers, dy, dx, dd, + friction_np, f_min, use_friction, + cellsize_x, cellsize_y) + path_data = cupy.asarray(path_img) - if start_py != NONE: - neighbor_ys, neighbor_xs = _neighborhood_structure(connectivity) - _a_star_search(surface.data, path_img, start_py, start_px, - goal_py, goal_px, barriers, neighbor_ys, neighbor_xs) - - path_agg = xr.DataArray(path_img, + else: + # numpy path (existing, unchanged) + if use_friction: + friction_data = np.asarray(friction.data, dtype=np.float64) + mask = np.isfinite(friction_data) & (friction_data > 0) + if not np.any(mask): + raise ValueError("friction has no positive finite values") + f_min = float(np.min(friction_data[mask])) + else: + friction_data = np.ones((h, w), dtype=np.float64) + f_min = 1.0 + + path_img = np.full(surface.shape, np.nan, dtype=np.float64) + if start_py != NONE: + _a_star_search(surface_data, path_img, start_py, start_px, + goal_py, goal_px, barriers, dy, dx, dd, + friction_data, f_min, use_friction, + cellsize_x, cellsize_y) + path_data = path_img + + path_agg = xr.DataArray(path_data, coords=surface.coords, dims=surface.dims, attrs=surface.attrs) diff --git a/xrspatial/tests/test_pathfinding.py b/xrspatial/tests/test_pathfinding.py index 0f718e12..d77e71fb 100644 --- a/xrspatial/tests/test_pathfinding.py +++ b/xrspatial/tests/test_pathfinding.py @@ -1,8 +1,18 @@ import numpy as np import pytest +try: + import dask.array as da +except ImportError: + da = None + from xrspatial import a_star_search -from xrspatial.tests.general_checks import create_test_raster, general_output_checks +from xrspatial.cost_distance import cost_distance +from xrspatial.tests.general_checks import ( + create_test_raster, general_output_checks, + cuda_and_cupy_available, dask_array_available, +) +from xrspatial.utils import has_cuda_and_cupy, has_dask_array @pytest.fixture @@ -33,20 +43,27 @@ def input_data_with_nans(): @pytest.fixture def result_8_connectivity(): + # cellsize=0.5: diagonal=sqrt(0.5)β‰ˆ0.7071, cardinal=0.5 + # Path: (0,2)β†’(1,1)β†’(2,1)β†’(3,1) + # Costs: 0, sqrt(0.5), sqrt(0.5)+0.5, sqrt(0.5)+1.0 + s = np.sqrt(0.5) expected_result = np.array([[np.nan, np.nan, 0., np.nan], - [np.nan, 1.41421356, np.nan, np.nan], - [np.nan, 2.41421356, np.nan, np.nan], - [np.nan, 3.41421356, np.nan, np.nan], + [np.nan, s, np.nan, np.nan], + [np.nan, s + 0.5, np.nan, np.nan], + [np.nan, s + 1.0, np.nan, np.nan], [np.nan, np.nan, np.nan, np.nan]]) return expected_result @pytest.fixture def result_4_connectivity(): - expected_result = np.array([[np.nan, 1, 0., np.nan], - [np.nan, 2, np.nan, np.nan], - [np.nan, 3, np.nan, np.nan], - [np.nan, 4, np.nan, np.nan], + # cellsize=0.5: cardinal=0.5 + # Path: (0,2)β†’(0,1)β†’(1,1)β†’(2,1)β†’(3,1) + # Costs: 0, 0.5, 1.0, 1.5, 2.0 + expected_result = np.array([[np.nan, 0.5, 0., np.nan], + [np.nan, 1.0, np.nan, np.nan], + [np.nan, 1.5, np.nan, np.nan], + [np.nan, 2.0, np.nan, np.nan], [np.nan, np.nan, np.nan, np.nan]]) return expected_result @@ -138,3 +155,413 @@ def test_a_star_search_connectivity( agg, start, goal, barriers, 'lon', 'lat', snap_start=True, snap_goal=True, connectivity=4 ) np.testing.assert_allclose(path_agg_4, result_4_connectivity, equal_nan=True) + + +# ----------------------------------------------------------------------- +# Helper for multi-backend test rasters +# ----------------------------------------------------------------------- + +def _make_raster(data, dims=None, res=None, backend='numpy', chunks=(3, 3)): + """Build a simple DataArray for weighted A* tests.""" + if dims is None: + dims = ['y', 'x'] + if res is None: + res = (1.0, 1.0) + import xarray as xr + h, w = data.shape + raster = xr.DataArray( + data.astype(np.float64), + dims=dims, + attrs={'res': res}, + ) + raster[dims[0]] = np.linspace((h - 1) * res[0], 0, h) + raster[dims[1]] = np.linspace(0, (w - 1) * res[1], w) + + if 'dask' in backend and da is not None: + raster.data = da.from_array(raster.data, chunks=chunks) + if 'cupy' in backend and has_cuda_and_cupy(): + import cupy + if isinstance(raster.data, da.Array): + raster.data = raster.data.map_blocks(cupy.asarray) + else: + raster.data = cupy.asarray(raster.data) + + return raster + + +# ----------------------------------------------------------------------- +# Weighted A* tests β€” parametrized for numpy and dask +# ----------------------------------------------------------------------- + +_backends = ['numpy'] +if has_dask_array(): + _backends.append('dask+numpy') + + +@pytest.mark.parametrize("backend", _backends) +def test_uniform_friction_matches_no_friction(backend): + """Uniform friction=1 should give the same path and costs as no friction.""" + data = np.ones((5, 5)) + agg = _make_raster(data, backend=backend) + friction_agg = _make_raster(np.ones((5, 5)), backend=backend) + + # coordinate space: y goes 4..0, x goes 0..4 + start = (4.0, 0.0) # pixel (0,0) + goal = (0.0, 4.0) # pixel (4,4) + + path_no_friction = a_star_search(agg, start, goal) + path_with_friction = a_star_search(agg, start, goal, friction=friction_agg) + + np.testing.assert_allclose( + np.asarray(path_with_friction.values), + np.asarray(path_no_friction.values), + equal_nan=True, atol=1e-10, + ) + + +@pytest.mark.parametrize("backend", _backends) +def test_high_friction_detour(backend): + """A* should route around a high-friction column.""" + # 3x5 grid, cellsize=1 + data = np.ones((3, 5)) + friction_data = np.ones((3, 5)) + friction_data[0, 2] = 100.0 # high friction in direct path + friction_data[1, 2] = 100.0 + + agg = _make_raster(data, backend=backend) + friction_agg = _make_raster(friction_data, backend=backend) + + # start=pixel(0,0), goal=pixel(0,4) + # coords: y=[2,1,0], x=[0,1,2,3,4] + start = (2.0, 0.0) + goal = (2.0, 4.0) + + # Without friction: direct straight path (0,0)β†’(0,1)β†’(0,2)β†’(0,3)β†’(0,4) + path_no = a_star_search(agg, start, goal) + assert np.isfinite(np.asarray(path_no.values)[0, 2]) # goes through column 2 + + # With friction: should detour around column 2 + path_fr = a_star_search(agg, start, goal, friction=friction_agg) + vals = np.asarray(path_fr.values) + assert np.isnan(vals[0, 2]) # does NOT go through (0,2) + # detour should go through row 2 (pixel row 2) + assert np.isfinite(vals[2, 2]) # goes through (2,2) instead + + # Weighted path should be cheaper than going through high-friction zone + assert np.nanmax(vals) < 10.0 + + +@pytest.mark.parametrize("backend", _backends) +def test_friction_barrier_nan(backend): + """NaN friction blocks passage, forcing a detour.""" + data = np.ones((3, 5)) + friction_data = np.ones((3, 5)) + friction_data[0, 2] = np.nan # impassable + + agg = _make_raster(data, backend=backend) + friction_agg = _make_raster(friction_data, backend=backend) + + start = (2.0, 0.0) # pixel (0,0) + goal = (2.0, 4.0) # pixel (0,4) + + path = a_star_search(agg, start, goal, friction=friction_agg) + vals = np.asarray(path.values) + # Path must not go through (0,2) + assert np.isnan(vals[0, 2]) + # But a path should still exist + assert np.isfinite(vals[0, 4]) + + +@pytest.mark.parametrize("backend", _backends) +def test_friction_barrier_zero(backend): + """Zero friction blocks passage, forcing a detour.""" + data = np.ones((3, 5)) + friction_data = np.ones((3, 5)) + friction_data[0, 2] = 0.0 # impassable + + agg = _make_raster(data, backend=backend) + friction_agg = _make_raster(friction_data, backend=backend) + + start = (2.0, 0.0) + goal = (2.0, 4.0) + + path = a_star_search(agg, start, goal, friction=friction_agg) + vals = np.asarray(path.values) + assert np.isnan(vals[0, 2]) + assert np.isfinite(vals[0, 4]) + + +@pytest.mark.parametrize("backend", _backends) +def test_a_star_matches_cost_distance(backend): + """A* path cost at goal should equal cost_distance value at that pixel.""" + # 5x5 grid with single source at (0,0), non-uniform friction + source = np.zeros((5, 5)) + source[0, 0] = 1.0 # source/start pixel + + friction_data = np.ones((5, 5)) + friction_data[1, 1] = 3.0 + friction_data[2, 2] = 2.0 + + surface = _make_raster(np.ones((5, 5)), backend=backend) + source_r = _make_raster(source, backend=backend) + friction_r = _make_raster(friction_data, backend=backend) + + # Run cost_distance from source at (0,0) β€” always numpy for comparison + source_np = _make_raster(source, backend='numpy') + friction_np = _make_raster(friction_data, backend='numpy') + cd = cost_distance(source_np, friction_np) + cd_val = cd.values + + # Run A* from (0,0) to (4,4) + # coords: y=[4,3,2,1,0], x=[0,1,2,3,4] + start = (4.0, 0.0) # pixel (0,0) + goal = (0.0, 4.0) # pixel (4,4) + + path = a_star_search(surface, start, goal, friction=friction_r) + a_star_cost_at_goal = float(np.asarray(path.values)[4, 4]) + + # cost_distance gives minimum cost from source to every pixel + cd_cost_at_goal = float(cd_val[4, 4]) + + np.testing.assert_allclose(a_star_cost_at_goal, cd_cost_at_goal, rtol=1e-5) + + +@pytest.mark.parametrize("backend", _backends) +def test_analytic_3x3_weighted(backend): + """Hand-computed costs on a 3x3 grid with known friction.""" + # Surface all 1s, friction has high center + data = np.ones((3, 3)) + friction_data = np.array([ + [1.0, 2.0, 1.0], + [1.0, 1.0, 1.0], + [1.0, 1.0, 1.0], + ]) + + agg = _make_raster(data, backend=backend) + friction_agg = _make_raster(friction_data, backend=backend) + + # start=pixel(0,0), goal=pixel(0,2) + # coords: y=[2,1,0], x=[0,1,2] + start = (2.0, 0.0) + goal = (2.0, 2.0) + + path = a_star_search(agg, start, goal, friction=friction_agg) + vals = np.asarray(path.values) + + # Direct path through (0,1) with friction=2: + # (0,0)β†’(0,1): 1*(1+2)/2 = 1.5 + # (0,1)β†’(0,2): 1*(2+1)/2 = 1.5 + # total = 3.0 + # + # Detour through (1,1) with friction=1: + # (0,0)β†’(1,1): sqrt(2)*(1+1)/2 = sqrt(2) β‰ˆ 1.4142 + # (1,1)β†’(0,2): sqrt(2)*(1+1)/2 = sqrt(2) β‰ˆ 1.4142 + # total = 2*sqrt(2) β‰ˆ 2.8284 + # + # A* picks the cheaper detour + expected_cost_at_goal = 2.0 * np.sqrt(2.0) + np.testing.assert_allclose(vals[0, 2], expected_cost_at_goal, rtol=1e-5) + + # Path should not go through (0,1) β€” the high-friction direct route + assert np.isnan(vals[0, 1]) + # Path should go through (1,1) β€” the low-friction detour + np.testing.assert_allclose(vals[1, 1], np.sqrt(2.0), rtol=1e-5) + + +# ----------------------------------------------------------------------- +# Dask-specific tests +# ----------------------------------------------------------------------- + +@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array") +@pytest.mark.parametrize("chunks", [(2, 2), (3, 3), (5, 5)]) +def test_dask_matches_numpy(chunks): + """Dask A* with various chunk sizes should match numpy results exactly.""" + data = np.ones((5, 5)) + friction_data = np.ones((5, 5)) + friction_data[1, 1] = 3.0 + friction_data[2, 2] = 2.0 + + agg_np = _make_raster(data, backend='numpy') + friction_np = _make_raster(friction_data, backend='numpy') + + agg_dask = _make_raster(data, backend='dask+numpy', chunks=chunks) + friction_dask = _make_raster(friction_data, backend='dask+numpy', + chunks=chunks) + + start = (4.0, 0.0) + goal = (0.0, 4.0) + + path_np = a_star_search(agg_np, start, goal, friction=friction_np) + path_dask = a_star_search(agg_dask, start, goal, friction=friction_dask) + + np.testing.assert_allclose( + np.asarray(path_dask.values), + path_np.values, + equal_nan=True, atol=1e-10, + ) + + +@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array") +def test_dask_no_large_numpy_arrays(): + """Dask path should not materialise full-size numpy arrays.""" + height, width = 50, 60 + data = np.ones((height, width)) + friction_data = np.ones((height, width)) + + agg = _make_raster(data, backend='dask+numpy', chunks=(25, 30)) + friction_agg = _make_raster(friction_data, backend='dask+numpy', + chunks=(25, 30)) + + start = (float(height - 1), 0.0) + goal = (0.0, float(width - 1)) + + # Track large numpy allocations + original_full = np.full + large_allocs = [] + + def tracking_full(shape, *args, **kwargs): + result = original_full(shape, *args, **kwargs) + if hasattr(shape, '__len__'): + total = 1 + for s in shape: + total *= s + else: + total = shape + if total >= height * width: + large_allocs.append(('full', shape)) + return result + + from unittest.mock import patch + with patch('numpy.full', side_effect=tracking_full): + result = a_star_search(agg, start, goal, friction=friction_agg) + + # Result should be dask-backed + assert isinstance(result.data, da.Array) + + # No full-size arrays should have been allocated + assert len(large_allocs) == 0, ( + f"Unexpected large allocations: {large_allocs}") + + +@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array") +def test_dask_snap_raises(): + """snap_start/snap_goal should raise ValueError for dask arrays.""" + data = np.ones((5, 5)) + agg = _make_raster(data, backend='dask+numpy', chunks=(3, 3)) + start = (4.0, 0.0) + goal = (0.0, 4.0) + + with pytest.raises(ValueError, match="snap_start is not supported"): + a_star_search(agg, start, goal, snap_start=True) + + with pytest.raises(ValueError, match="snap_goal is not supported"): + a_star_search(agg, start, goal, snap_goal=True) + + +@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array") +def test_dask_no_friction(): + """Dask A* without friction should match numpy.""" + data = np.ones((5, 5)) + agg_np = _make_raster(data, backend='numpy') + agg_dask = _make_raster(data, backend='dask+numpy', chunks=(3, 3)) + + start = (4.0, 0.0) + goal = (0.0, 4.0) + + path_np = a_star_search(agg_np, start, goal) + path_dask = a_star_search(agg_dask, start, goal) + + np.testing.assert_allclose( + np.asarray(path_dask.values), + path_np.values, + equal_nan=True, atol=1e-10, + ) + + +# ----------------------------------------------------------------------- +# CuPy tests +# ----------------------------------------------------------------------- + +@cuda_and_cupy_available +def test_cupy_matches_numpy(): + """CuPy path should produce identical results to numpy.""" + data = np.ones((5, 5)) + friction_data = np.ones((5, 5)) + friction_data[1, 1] = 3.0 + friction_data[2, 2] = 2.0 + + agg_np = _make_raster(data, backend='numpy') + friction_np = _make_raster(friction_data, backend='numpy') + + agg_cupy = _make_raster(data, backend='cupy') + friction_cupy = _make_raster(friction_data, backend='cupy') + + start = (4.0, 0.0) + goal = (0.0, 4.0) + + path_np = a_star_search(agg_np, start, goal, friction=friction_np) + path_cupy = a_star_search(agg_cupy, start, goal, friction=friction_cupy) + + np.testing.assert_allclose( + path_cupy.data.get(), + path_np.values, + equal_nan=True, atol=1e-10, + ) + + +@cuda_and_cupy_available +def test_cupy_no_friction(): + """CuPy without friction should match numpy.""" + data = np.ones((5, 5)) + agg_np = _make_raster(data, backend='numpy') + agg_cupy = _make_raster(data, backend='cupy') + + start = (4.0, 0.0) + goal = (0.0, 4.0) + + path_np = a_star_search(agg_np, start, goal) + path_cupy = a_star_search(agg_cupy, start, goal) + + np.testing.assert_allclose( + path_cupy.data.get(), + path_np.values, + equal_nan=True, atol=1e-10, + ) + + +# ----------------------------------------------------------------------- +# Dask + CuPy tests +# ----------------------------------------------------------------------- + +@cuda_and_cupy_available +@pytest.mark.skipif(not has_dask_array(), reason="Requires dask.Array") +def test_dask_cupy_matches_numpy(): + """Dask+CuPy path should produce identical results to numpy.""" + data = np.ones((5, 5)) + friction_data = np.ones((5, 5)) + friction_data[1, 1] = 3.0 + friction_data[2, 2] = 2.0 + + agg_np = _make_raster(data, backend='numpy') + friction_np = _make_raster(friction_data, backend='numpy') + + agg_dc = _make_raster(data, backend='dask+cupy', chunks=(3, 3)) + friction_dc = _make_raster(friction_data, backend='dask+cupy', + chunks=(3, 3)) + + start = (4.0, 0.0) + goal = (0.0, 4.0) + + path_np = a_star_search(agg_np, start, goal, friction=friction_np) + path_dc = a_star_search(agg_dc, start, goal, friction=friction_dc) + + # dask+cupy: compute dask graph, then .get() the cupy array to numpy + dc_computed = path_dc.data.compute() + if hasattr(dc_computed, 'get'): + dc_computed = dc_computed.get() + + np.testing.assert_allclose( + dc_computed, + path_np.values, + equal_nan=True, atol=1e-10, + )