"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "blue | green |
red | nir |
\n",
- "
"
- ],
- "text/plain": [
- "<xarray.DataArray 'true_color' (y: 1024, x: 1024, band: 4)> Size: 4MB\n", - "array([[[56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " ...,\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0]],\n", - "\n", - " [[56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " ...,\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0]],\n", - "\n", - " [[56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " ...,\n", - "...\n", - " ...,\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0]],\n", - "\n", - " [[56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " ...,\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0]],\n", - "\n", - " [[56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " ...,\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0],\n", - " [56, 56, 56, 0]]], shape=(1024, 1024, 4), dtype=uint8)\n", - "Coordinates:\n", - " * y (y) float64 8kB 6.94e+06 6.94e+06 6.94e+06 ... 6.698e+06 6.697e+06\n", - " * x (x) float64 8kB 5.356e+05 5.358e+05 ... 7.758e+05 7.76e+05\n", - " * band (band) int64 32B 0 1 2 3\n", - "Attributes:\n", - " res: 30.0\n", - " x_range: (np.float64(535485.0), np.float64(776115.0))\n", - " y_range: (np.float64(6697185.0), np.float64(6940515.0))
nir | red | ndvi |
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- ""
- ],
- "text/plain": [
- "
"
- ],
- "text/plain": [
- "\n",
- "\n",
- "\n",
- " \n",
- "
\n",
- "
"
- ],
- "text/plain": [
- " zone mean max min sum std \\\n",
- "0 0 888.487944 4000.000000 0.000000 4.264715e+08 872.385158 \n",
- "1 11 1443.760982 1452.024715 1435.497249 2.887522e+03 8.263733 \n",
- "2 70 1435.677678 1435.677678 1435.677678 1.435678e+03 0.000000 \n",
- "\n",
- " var count \n",
- "0 761055.863901 479997.0 \n",
- "1 68.289285 2.0 \n",
- "2 0.000000 1.0 "
- ]
- },
- "execution_count": 5,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
+ "outputs": [],
"source": [
"from xrspatial import zonal_stats\n",
"\n",
"zones_agg.values = np.nan_to_num(zones_agg.values, copy=False).astype(int)\n",
- "zonal_stats(zones_agg, terrain)"
+ "zonal_stats(zones_agg, terrain)\n"
]
},
{
@@ -398,7 +280,7 @@
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": null,
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-23T01:08:24.736186Z",
@@ -407,72 +289,7 @@
"shell.execute_reply": "2026-01-23T01:08:24.774232Z"
}
},
- "outputs": [
- {
- "data": {
- "text/html": [
- "| \n", - " | zone | \n", - "mean | \n", - "max | \n", - "min | \n", - "sum | \n", - "std | \n", - "var | \n", - "count | \n", - "
|---|---|---|---|---|---|---|---|---|
| 0 | \n", - "0 | \n", - "888.487944 | \n", - "4000.000000 | \n", - "0.000000 | \n", - "4.264715e+08 | \n", - "872.385158 | \n", - "761055.863901 | \n", - "479997.0 | \n", - "
| 1 | \n", - "11 | \n", - "1443.760982 | \n", - "1452.024715 | \n", - "1435.497249 | \n", - "2.887522e+03 | \n", - "8.263733 | \n", - "68.289285 | \n", - "2.0 | \n", - "
| 2 | \n", - "70 | \n", - "1435.677678 | \n", - "1435.677678 | \n", - "1435.677678 | \n", - "1.435678e+03 | \n", - "0.000000 | \n", - "0.000000 | \n", - "1.0 | \n", - "
\n",
- "\n",
- "\n",
- " \n",
- "
\n",
- "
"
- ],
- "text/plain": [
- " zone elevation_change elevation_min elevation_max\n",
- "0 0 4000.000000 0.000000 4000.000000\n",
- "1 11 16.527466 1435.497249 1452.024715\n",
- "2 70 0.000000 1435.677678 1435.677678"
- ]
- },
- "execution_count": 6,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
+ "outputs": [],
"source": [
"custom_stats = dict(\n",
" elevation_change=lambda zone: zone.max() - zone.min(),\n",
@@ -480,14 +297,14 @@
" elevation_max=np.max,\n",
")\n",
"\n",
- "zonal_stats(zones=zones_agg, values=terrain, stats_funcs=custom_stats)"
+ "zonal_stats(zones=zones_agg, values=terrain, stats_funcs=custom_stats)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Here the zones are defined by line segments, but they can be any spatial pattern or, more specifically, any region computable as a Datashader aggregate."
+ "Here the zones are defined by line segments, but they can be any spatial pattern or, more specifically, any region computable as a Raster aggregate."
]
}
],
diff --git a/examples/Pathfinding_Austin_Road_Network.ipynb b/examples/Pathfinding_Austin_Road_Network.ipynb
index 0e5e565b3..9fdab101e 100644
--- a/examples/Pathfinding_Austin_Road_Network.ipynb
+++ b/examples/Pathfinding_Austin_Road_Network.ipynb
@@ -18,7 +18,7 @@
"\n",
"First, we'll need to import some packages: these include the basic array manipulation ones, \n",
"as well as some geospatial-focused ones.\n",
- "We'll also grab a few datashader functions for easy rendering."
+ "We'll also grab a few matplotlib functions for easy rendering."
]
},
{
@@ -34,15 +34,68 @@
},
"outputs": [],
"source": [
- "import datashader as ds\n",
- "from datashader.transfer_functions import shade, stack, dynspread, set_background\n",
"import geopandas\n",
"import numpy as np\n",
"import pandas as pd\n",
- "import spatialpandas\n",
"import xarray as xa\n",
"\n",
- "from xrspatial import a_star_search"
+ "from xrspatial import a_star_search\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "execution_count": null,
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib.colors import LinearSegmentedColormap, to_rgba\n",
+ "import numpy as np\n",
+ "\n",
+ "def _shade(arr, cmap=None, alpha=None, min_alpha=None, span=None, how='linear'):\n",
+ " arr = np.asarray(arr, dtype=np.float64)\n",
+ " finite = np.isfinite(arr)\n",
+ " if span is not None:\n",
+ " lo, hi = span\n",
+ " elif finite.any():\n",
+ " lo = np.nanmin(arr); hi = np.nanmax(arr)\n",
+ " else:\n",
+ " lo = hi = 0.0\n",
+ " norm = np.zeros_like(arr, dtype=np.float32)\n",
+ " if hi > lo:\n",
+ " norm = np.where(finite, (arr - lo) / (hi - lo), 0.0).astype(np.float32)\n",
+ " if isinstance(cmap, (list, tuple)):\n",
+ " cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)\n",
+ " elif cmap is None:\n",
+ " cmap = plt.get_cmap('terrain')\n",
+ " rgba = cmap(norm, bytes=True)\n",
+ " a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)\n",
+ " if min_alpha is not None:\n",
+ " a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)\n",
+ " rgba[..., 3] = a\n",
+ " return rgba\n",
+ "\n",
+ "def _stack(*layers):\n",
+ " base = layers[0].astype(np.float64)\n",
+ " for top in layers[1:]:\n",
+ " a = (top[..., 3:4] / 255.0)\n",
+ " rgb = top[..., :3] * a + base[..., :3] * (1 - a)\n",
+ " out_a = np.maximum(top[..., 3], base[..., 3])\n",
+ " base = np.concatenate([rgb, out_a[..., None]], axis=-1)\n",
+ " return np.clip(base, 0, 255).astype(np.uint8)\n",
+ "\n",
+ "def _show(rgba, bg=None, figsize=(10, 8)):\n",
+ " if bg is not None:\n",
+ " h, w = rgba.shape[:2]\n",
+ " bgc = to_rgba(bg)\n",
+ " bg_rgba = np.zeros((h, w, 4), dtype=np.uint8)\n",
+ " bg_rgba[..., :3] = (np.array(bgc[:3]) * 255).astype(np.uint8)\n",
+ " bg_rgba[..., 3] = 255\n",
+ " rgba = _stack(bg_rgba, rgba)\n",
+ " plt.figure(figsize=figsize)\n",
+ " plt.imshow(rgba)\n",
+ " plt.axis('off')\n",
+ " plt.show()\n"
]
},
{
@@ -61,7 +114,7 @@
"metadata": {},
"source": [
"We'll start by opening the shapefile, transforming the crs (coordinate reference system) to the commonly-used longitude/latitude, \n",
- "and, after a quick clean-up, transforming it to a spatialpandas GeoDataFrame.\n",
+ "and.\n",
"\n",
"Now our data is ready to be aggregated to an xarray DataArray raster."
]
@@ -81,8 +134,7 @@
"source": [
"streets = geopandas.read_file(\"../xrspatial-examples/data/geo_export_9c395dda-0b29-41ec-89b4-a51a898f7104.shp\")\n",
"streets = streets.to_crs(\"EPSG:4326\")\n",
- "streets = streets.explode(\"geometry\").reset_index(drop=True)\n",
- "streets_spd = spatialpandas.GeoDataFrame(streets, geometry=\"geometry\")"
+ "streets = streets.explode(\"geometry\").reset_index(drop=True)\n"
]
},
{
@@ -93,7 +145,7 @@
"\n",
"To finish off our set-up:\n",
"- We'll define a study area, with xmin, xmax, ymin, and ymax; this set the x, y coordinates we'll be using in our aggregate.\n",
- "- We'll set up a datashader Canvas object, which provides an easy frame for setting up a new raster and aggregating data to it.\n",
+ "- We'll set up a matplotlib Canvas object, which provides an easy frame for setting up a new raster and aggregating data to it.\n",
"- Finally, we'll aggregate the streets data into a lines raster with Canvas.line.\n",
"\n",
"- We also set up the start and goal point (y, x) coordinates, and set up a DataFrame and aggregation for visualization.\n",
@@ -122,7 +174,7 @@
")\n",
"xrange = (xmin, xmax)\n",
"yrange = (ymin, ymax)\n",
- "xrange, yrange"
+ "xrange, yrange\n"
]
},
{
@@ -139,24 +191,32 @@
"outputs": [],
"source": [
"H, W = 600, 800\n",
- "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=xrange, y_range=yrange)\n",
+ "grid_xs = np.linspace(xrange[0], xrange[1], W)\n",
+ "grid_ys = np.linspace(yrange[1], xrange[0], H)\n",
+ "template = xa.DataArray(np.full((H, W), np.nan),\n",
+ " coords={'y': grid_ys, 'x': grid_xs}, dims=['y', 'x'])\n",
"\n",
- "street_agg = cvs.line(streets_spd, geometry=\"geometry\")\n",
- "street_shaded = dynspread(shade(street_agg, cmap=[\"salmon\"]))\n",
+ "# Rasterize the street lines onto the grid with the .xrs accessor\n",
+ "# (replaces ds.Canvas().line(streets_spd, geometry='geometry')).\n",
+ "from shapely.geometry import Point\n",
+ "street_agg = template.xrs.rasterize(streets, merge='last')\n",
+ "street_shaded = _shade(street_agg, cmap=['salmon', 'salmon'], min_alpha=255)\n",
"\n",
- "# Pick two locations\n",
+ "# Pick two locations: (lat, lon) pairs\n",
"start = (30.08214069, -97.73662282)\n",
"goal = (30.17656606, -97.63753489)\n",
"\n",
- "start_df = pd.DataFrame({\"x\": [start[1]], \"y\": [start[0]]})\n",
- "start_agg = cvs.points(start_df, \"x\", \"y\")\n",
- "start_shaded = dynspread(shade(start_agg, cmap=[\"red\"]), threshold=1, max_px=5)\n",
+ "start_gdf = geopandas.GeoDataFrame(\n",
+ " {'id': [1]}, geometry=[Point(start[1], start[0])], crs='EPSG:4326')\n",
+ "start_agg = template.xrs.rasterize(start_gdf, column='id', merge='last')\n",
+ "start_shaded = _shade(start_agg, cmap=['red', 'red'], min_alpha=255)\n",
"\n",
- "goal_df = pd.DataFrame({\"x\": [goal[1]], \"y\": [goal[0]]})\n",
- "goal_agg = cvs.points(goal_df, \"x\", \"y\")\n",
- "goal_shaded = dynspread(shade(goal_agg, cmap=[\"lime\"]), threshold=1, max_px=5)\n",
+ "goal_gdf = geopandas.GeoDataFrame(\n",
+ " {'id': [1]}, geometry=[Point(goal[1], goal[0])], crs='EPSG:4326')\n",
+ "goal_agg = template.xrs.rasterize(goal_gdf, column='id', merge='last')\n",
+ "goal_shaded = _shade(goal_agg, cmap=['lime', 'lime'], min_alpha=255)\n",
"\n",
- "set_background(stack(street_shaded, start_shaded, goal_shaded), \"black\")"
+ "_show(_stack(street_shaded, start_shaded, goal_shaded), bg='black')\n"
]
},
{
@@ -196,8 +256,8 @@
" street_agg, start, goal, barriers=[0], snap_start=True, snap_goal=True\n",
")\n",
"\n",
- "path_shaded = dynspread(shade(path_agg, cmap=[\"green\"]), threshold=1, max_px=1)\n",
- "set_background(stack(street_shaded, path_shaded, start_shaded, goal_shaded), \"black\")"
+ "path_shaded = _shade(path_agg, cmap=['green', 'green'], min_alpha=255)\n",
+ "_show(_stack(street_shaded, path_shaded, start_shaded, goal_shaded), bg='black')\n"
]
},
{
diff --git a/examples/animated_hillshade.py b/examples/animated_hillshade.py
index 9298cc017..b12036306 100644
--- a/examples/animated_hillshade.py
+++ b/examples/animated_hillshade.py
@@ -1,20 +1,34 @@
+"""Animated hillshade example -- matplotlib / PIL rendering (no datashader).
+
+Generates a synthetic terrain, places trees and water on it, and renders a
+rotating-hillshade GIF plus a composited build-up GIF.
+
+Requires: numpy, xarray, matplotlib, pillow, xrspatial.
+"""
from functools import partial
-import datashader as ds
import numpy as np
-from datashader.colors import Elevation
-from datashader.transfer_functions import shade, stack
+import xarray as xr
+import matplotlib.cm as cm
+from matplotlib.colors import LinearSegmentedColormap
+from PIL import Image
from xrspatial import bump, generate_terrain, hillshade, mean
W = 600
H = 400
-cvs = ds.Canvas(plot_width=W, plot_height=H,
- x_range=(-20e6, 20e6),
- y_range=(-20e6, 20e6))
+x_range = (-20e6, 20e6)
+y_range = (-20e6, 20e6)
+
+xs = np.linspace(x_range[0], x_range[1], W)
+ys = np.linspace(y_range[1], y_range[0], H) # origin upper-left to match rasters
+template = xr.DataArray(
+ np.zeros((H, W), dtype=np.float32),
+ coords={'y': ys, 'x': xs}, dims=['y', 'x'], name='terrain',
+)
-terrain = generate_terrain(canvas=cvs)
+terrain = generate_terrain(template, x_range=x_range, y_range=y_range)
def heights(locations, src, src_range, height=20):
@@ -34,10 +48,10 @@ def heights(locations, src, src_range, height=20):
src = terrain.data
trees = bump(W, H, count=T, height_func=partial(heights, src=src,
src_range=(1000, 1300), height=5))
-trees += bump(W, H, count=T//2, height_func=partial(
- heights, src=src, src_range=(1300, 1700), height=20))
-trees += bump(W, H, count=T//3, height_func=partial(
- heights, src=src, src_range=(1700, 2000), height=5))
+trees += bump(W, H, count=T // 2, height_func=partial(
+ heights, src=src, src_range=(1300, 1700), height=20))
+trees += bump(W, H, count=T // 3, height_func=partial(
+ heights, src=src, src_range=(1700, 2000), height=5))
tree_colorize = trees.copy()
tree_colorize.data[tree_colorize.data == 0] = np.nan
@@ -50,84 +64,134 @@ def heights(locations, src, src_range, height=20):
water.data[water.data == LAND_CONSTANT] = np.nan
-def create_map(azimuth):
-
- global cvs
- global terrain
- global water
- global trees
+# Elevation-style colormap (replaces datashader.colors.Elevation).
+_ELEV_STOPS = [
+ (0.00, (0.00, 0.30, 0.60)), # deep water
+ (0.20, (0.20, 0.60, 0.85)), # shallow water
+ (0.30, (0.55, 0.78, 0.92)), # shoreline
+ (0.32, (0.40, 0.73, 0.45)), # beach / lowland
+ (0.45, (0.30, 0.55, 0.25)), # forested
+ (0.60, (0.50, 0.50, 0.30)), # shrubland
+ (0.75, (0.65, 0.55, 0.40)), # highland
+ (0.90, (0.80, 0.80, 0.78)), # bare rock
+ (1.00, (1.00, 1.00, 1.00)), # snow
+]
+ELEVATION = LinearSegmentedColormap.from_list('elevation',
+ [s[1] for s in _ELEV_STOPS],
+ N=256)
+
+
+def _normalize(arr):
+ """Scale a finite-valued array to [0, 1], NaN -> 0 (transparent later)."""
+ finite = np.isfinite(arr)
+ if not finite.any():
+ return np.zeros_like(arr, dtype=np.float32)
+ lo = np.nanmin(arr)
+ hi = np.nanmax(arr)
+ if hi <= lo:
+ return np.zeros_like(arr, dtype=np.float32)
+ out = np.where(finite, (arr - lo) / (hi - lo), 0.0)
+ return out.astype(np.float32)
+
+
+def _rgba(arr, cmap, alpha=None):
+ """Map a 2D array through a matplotlib colormap to an (H, W, 4) uint8 image.
+
+ NaN pixels become fully transparent (alpha=0). If `alpha` (0-255 int) is
+ given, every finite pixel gets that alpha instead of 255.
+ """
+ norm = _normalize(arr)
+ rgba = cm.ScalarMappable(cmap=cmap).to_rgba(norm, bytes=True)
+ # to_rgba ignores nan by setting alpha 0; restore finite alpha if requested
+ if alpha is not None:
+ finite = np.isfinite(arr)
+ rgba[..., 3] = np.where(finite, alpha, 0).astype(np.uint8)
+ else:
+ finite = np.isfinite(arr)
+ rgba[..., 3] = np.where(finite, 255, 0).astype(np.uint8)
+ return rgba
+
+
+def _composite(layers):
+ """Alpha-composite a list of (H, W, 4) uint8 arrays (bottom-to-top)."""
+ base = layers[0].astype(np.float64)
+ for top in layers[1:]:
+ a = (top[..., 3:4] / 255.0)
+ rgb = top[..., :3] * a + base[..., :3] * (1 - a)
+ out_a = np.maximum(top[..., 3:4], base[..., 3:4])
+ base = np.concatenate([rgb, out_a], axis=-1)
+ return np.clip(base, 0, 255).astype(np.uint8)
+
+
+def _to_pil(rgba):
+ return Image.fromarray(rgba, 'RGBA').convert('RGB')
- img = stack(shade(terrain, cmap=Elevation, how='linear'),
- shade(water, cmap=['aqua', 'white']),
- shade(hillshade(terrain + trees, azimuth=azimuth),
- cmap=['black', 'white'], how='linear', alpha=128),
- shade(tree_colorize, cmap='limegreen', how='linear')
- )
+def create_map(azimuth):
+ water_cmap = LinearSegmentedColormap.from_list('water',
+ ['aqua', 'white'], N=256)
+ layers = [
+ _rgba(terrain, ELEVATION, alpha=255),
+ _rgba(water, water_cmap, alpha=255),
+ _rgba(hillshade(terrain + trees, azimuth=azimuth),
+ LinearSegmentedColormap.from_list('hs', ['black', 'white'], N=256),
+ alpha=128),
+ _rgba(tree_colorize, LinearSegmentedColormap.from_list(
+ 'trees', ['limegreen', 'limegreen'], N=256), alpha=255),
+ ]
+ img = _to_pil(_composite(layers))
print('image created')
-
- return img.to_pil()
+ return img
def create_map2():
-
- global cvs
- global terrain
- global water
- global trees
-
- img = stack(shade(terrain, cmap=['black', 'white'], how='linear'))
-
- yield img.to_pil()
-
- img = stack(shade(terrain, cmap=Elevation, how='linear'))
-
- yield img.to_pil()
-
- img = stack(shade(terrain, cmap=Elevation, how='linear'),
- shade(hillshade(terrain, azimuth=210),
- cmap=['black', 'white'], how='linear', alpha=128),
- )
-
- yield img.to_pil()
-
- img = stack(shade(terrain, cmap=Elevation, how='linear'),
- shade(water, cmap=['aqua', 'white']),
- shade(hillshade(terrain, azimuth=210),
- cmap=['black', 'white'], how='linear', alpha=128),
- )
-
- yield img.to_pil()
-
- img = stack(shade(terrain, cmap=Elevation, how='linear'),
- shade(water, cmap=['aqua', 'white']),
- shade(hillshade(terrain + trees, azimuth=210),
- cmap=['black', 'white'], how='linear', alpha=128),
- shade(tree_colorize, cmap='limegreen', how='linear')
- )
-
- yield img.to_pil()
- yield img.to_pil()
- yield img.to_pil()
- yield img.to_pil()
+ water_cmap = LinearSegmentedColormap.from_list('water',
+ ['aqua', 'white'], N=256)
+ bw = LinearSegmentedColormap.from_list('hs', ['black', 'white'], N=256)
+
+ img = _to_pil(_rgba(terrain, bw, alpha=255))
+ yield img
+
+ img = _to_pil(_rgba(terrain, ELEVATION, alpha=255))
+ yield img
+
+ img = _to_pil(_composite([
+ _rgba(terrain, ELEVATION, alpha=255),
+ _rgba(hillshade(terrain, azimuth=210), bw, alpha=128),
+ ]))
+ yield img
+
+ img = _to_pil(_composite([
+ _rgba(terrain, ELEVATION, alpha=255),
+ _rgba(water, water_cmap, alpha=255),
+ _rgba(hillshade(terrain, azimuth=210), bw, alpha=128),
+ ]))
+ yield img
+
+ img = _to_pil(_composite([
+ _rgba(terrain, ELEVATION, alpha=255),
+ _rgba(water, water_cmap, alpha=255),
+ _rgba(hillshade(terrain + trees, azimuth=210), bw, alpha=128),
+ _rgba(tree_colorize, LinearSegmentedColormap.from_list(
+ 'trees', ['limegreen', 'limegreen'], N=256), alpha=255),
+ ]))
+ yield img
+ yield img
+ yield img
+ yield img
def gif1():
-
images = []
-
for i in np.linspace(0, 360, 6):
images.append(create_map(int(i)))
-
images[0].save('animated_hillshade.gif',
save_all=True, append_images=images[1:],
optimize=False, duration=5000, loop=0)
def gif2():
-
images = list(create_map2())
-
images[0].save('composite_map.gif',
save_all=True, append_images=images[1:],
optimize=False, duration=1000, loop=0)
diff --git a/examples/pharmacy-deserts.ipynb b/examples/pharmacy-deserts.ipynb
index 88e0c45fd..daf5ed82e 100644
--- a/examples/pharmacy-deserts.ipynb
+++ b/examples/pharmacy-deserts.ipynb
@@ -15,7 +15,7 @@
"source": [
"### Setup:\n",
"\n",
- "First, we'll import some basic data and geospatial data-related packages. We'll also bring in some datashader functions for aggregation and rendering. And, finally, we'll import the Xarray-spatial functions we'll be using. "
+ "First, we'll import some basic data and geospatial data-related packages. We'll also bring in some matplotlib functions for aggregation and rendering. And, finally, we'll import the Xarray-spatial functions we'll be using. "
]
},
{
@@ -30,7 +30,72 @@
}
},
"outputs": [],
- "source": "import datashader as ds\nfrom datashader.colors import inferno\nfrom datashader.transfer_functions import set_background\nfrom datashader.transfer_functions import shade\nfrom datashader.transfer_functions import stack\nimport geopandas as gpd\nimport numpy as np\nimport pandas as pd\nfrom spatialpandas import GeoDataFrame\n\nimport xrspatial"
+ "source": [
+ "import geopandas as gpd\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "\n",
+ "import xrspatial\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "metadata": {},
+ "execution_count": null,
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib.colors import LinearSegmentedColormap, to_rgba\n",
+ "import numpy as np\n",
+ "\n",
+ "def _shade(arr, cmap=None, alpha=None, min_alpha=None, span=None, how='linear'):\n",
+ " arr = np.asarray(arr, dtype=np.float64)\n",
+ " finite = np.isfinite(arr)\n",
+ " if span is not None:\n",
+ " lo, hi = span\n",
+ " elif finite.any():\n",
+ " lo = np.nanmin(arr); hi = np.nanmax(arr)\n",
+ " else:\n",
+ " lo = hi = 0.0\n",
+ " norm = np.zeros_like(arr, dtype=np.float32)\n",
+ " if hi > lo:\n",
+ " norm = np.where(finite, (arr - lo) / (hi - lo), 0.0).astype(np.float32)\n",
+ " if isinstance(cmap, (list, tuple)):\n",
+ " cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)\n",
+ " elif cmap is None:\n",
+ " cmap = plt.get_cmap('terrain')\n",
+ " rgba = cmap(norm, bytes=True)\n",
+ " a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)\n",
+ " if min_alpha is not None:\n",
+ " a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)\n",
+ " rgba[..., 3] = a\n",
+ " return rgba\n",
+ "\n",
+ "def _stack(*layers):\n",
+ " base = layers[0].astype(np.float64)\n",
+ " for top in layers[1:]:\n",
+ " a = (top[..., 3:4] / 255.0)\n",
+ " rgb = top[..., :3] * a + base[..., :3] * (1 - a)\n",
+ " out_a = np.maximum(top[..., 3], base[..., 3])\n",
+ " base = np.concatenate([rgb, out_a[..., None]], axis=-1)\n",
+ " return np.clip(base, 0, 255).astype(np.uint8)\n",
+ "\n",
+ "def _show(rgba, bg=None, figsize=(12, 6)):\n",
+ " if bg is not None:\n",
+ " h, w = rgba.shape[:2]\n",
+ " bgc = to_rgba(bg)\n",
+ " bg_rgba = np.zeros((h, w, 4), dtype=np.uint8)\n",
+ " bg_rgba[..., :3] = (np.array(bgc[:3]) * 255).astype(np.uint8)\n",
+ " bg_rgba[..., 3] = 255\n",
+ " rgba = _stack(bg_rgba, rgba)\n",
+ " plt.figure(figsize=figsize)\n",
+ " plt.imshow(rgba)\n",
+ " plt.axis('off')\n",
+ " plt.show()\n",
+ "\n",
+ "inferno = plt.get_cmap('inferno')\n",
+ "Set1 = plt.get_cmap('Set1')\n"
+ ]
},
{
"cell_type": "markdown",
@@ -68,7 +133,7 @@
"coords = pharmacy_df[\"CalcLocation\"].str.split(\",\", expand=True)\n",
"pharmacy_df[\"y\"] = np.array(coords[0], dtype=\"float64\")\n",
"pharmacy_df[\"x\"] = np.array(coords[1], dtype=\"float64\")\n",
- "pharmacy_df"
+ "pharmacy_df\n"
]
},
{
@@ -94,12 +159,11 @@
"source": [
"# Load Census Block Groups and Calculate Percent over 65 years-old\n",
"blockgroup_df = gpd.read_file(\"../xrspatial-examples/data/USA_Block_Groups.shp\")\n",
- "blockgroup_df = GeoDataFrame(blockgroup_df, geometry=\"geometry\")\n",
"blockgroup_df[\"ABOVE_65\"] = blockgroup_df[[\"AGE_65_74\", \"AGE_75_84\", \"AGE_85_UP\"]].sum(\n",
" axis=1\n",
")\n",
"blockgroup_df[\"PCT_ABOVE_65\"] = blockgroup_df[\"ABOVE_65\"] / blockgroup_df[\"POP2010\"]\n",
- "blockgroup_df"
+ "blockgroup_df\n"
]
},
{
@@ -124,9 +188,8 @@
"source": [
"# Load Census County and Calculate Percent over 65 years-old\n",
"county_df = gpd.read_file(\"../xrspatial-examples/data/USA_Counties.shp\")\n",
- "county_df = GeoDataFrame(county_df, geometry=\"geometry\")\n",
"county_df[\"ZONE_ID\"] = county_df[\"OBJECTID\"].astype(np.int16)\n",
- "county_df"
+ "county_df\n"
]
},
{
@@ -145,10 +208,10 @@
"The first raster we'll set up is the counties mask; we'll use this one to make our other rasters fit with county shapes.\n",
"\n",
"- We set a range for for our x and y coordinates in longitude and latitude values.\n",
- "- We use datashader's Canvas object as an easy-to-use frame to build up a raster with and aggregate data into.\n",
+ "- We use matplotlib's Canvas object as an easy-to-use frame to build up a raster with and aggregate data into.\n",
"- Finally, we aggregate the county shapes - polygons - into that raster. \n",
"\n",
- "A quick shade and transformation with datashader functions visualizes our mask set-up."
+ "A quick shade and transformation with matplotlib functions visualizes our mask set-up."
]
},
{
@@ -170,10 +233,17 @@
"W = 1600\n",
"H = 800\n",
"\n",
- "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=x_range, y_range=y_range)\n",
+ "import xarray as xr\n",
+ "from shapely.geometry import Point\n",
+ "grid_xs = np.linspace(x_range[0], x_range[1], W)\n",
+ "grid_ys = np.linspace(y_range[1], y_range[0], H)\n",
+ "template = xr.DataArray(np.full((H, W), np.nan),\n",
+ " coords={'y': grid_ys, 'x': grid_xs}, dims=['y', 'x'])\n",
"\n",
- "county_mask = cvs.polygons(county_df, geometry=\"geometry\")\n",
- "set_background(shade(county_mask, cmap=\"#333333\", alpha=255), \"black\")"
+ "# Rasterize county polygons (no attribute) with the .xrs accessor\n",
+ "# (replaces ds.Canvas().polygons(county_df, geometry='geometry')).\n",
+ "county_mask = template.xrs.rasterize(county_df, merge='last')\n",
+ "_show(_shade(county_mask, cmap='#333333', alpha=255), bg='black')\n"
]
},
{
@@ -206,7 +276,25 @@
}
},
"outputs": [],
- "source": "pharmacy_raster = cvs.points(pharmacy_df, \"x\", \"y\")\nproximity_raster = pharmacy_raster.xrs.proximity(distance_metric=\"GREAT_CIRCLE\").where(\n county_mask\n)\nproximity_raster.data[~np.isfinite(proximity_raster.data)] = 0.0\n\nproximity_classifed = proximity_raster.xrs.natural_breaks(num_sample=20000, k=4).where(county_mask)\n\nimage_pharmacy = shade(proximity_classifed, cmap=inferno, alpha=255)\nimage_pharmacy = set_background(image_pharmacy, \"black\")\nimage_pharmacy"
+ "source": [
+ "# Build a GeoDataFrame of pharmacy Points (the CSV has x/y columns, no\n",
+ "# geometry) and rasterize with the .xrs accessor.\n",
+ "pharmacy_gdf = gpd.GeoDataFrame(\n",
+ " pharmacy_df[['CalcLocation']].copy(),\n",
+ " geometry=[Point(x, y) for x, y in zip(pharmacy_df['x'], pharmacy_df['y'])],\n",
+ " crs='EPSG:4326',\n",
+ ")\n",
+ "pharmacy_raster = template.xrs.rasterize(pharmacy_gdf, merge='last')\n",
+ "proximity_raster = pharmacy_raster.xrs.proximity(distance_metric=\"GREAT_CIRCLE\").where(\n",
+ " county_mask\n",
+ ")\n",
+ "proximity_raster.data[~np.isfinite(proximity_raster.data)] = 0.0\n",
+ "\n",
+ "proximity_classifed = proximity_raster.xrs.natural_breaks(num_sample=20000, k=4).where(county_mask)\n",
+ "\n",
+ "image_pharmacy = _shade(proximity_classifed, cmap=inferno, alpha=255)\n",
+ "_show(image_pharmacy, bg='black')\n"
+ ]
},
{
"cell_type": "markdown",
@@ -230,7 +318,20 @@
}
},
"outputs": [],
- "source": "age_raster = cvs.polygons(\n blockgroup_df, geometry=\"geometry\", agg=ds.mean(\"PCT_ABOVE_65\")\n)\nage_raster.data[~np.isfinite(age_raster.data)] = 0.0\n\nage_classifed = age_raster.xrs.natural_breaks(num_sample=20000, k=4).where(county_mask)\n\nage_image = shade(age_classifed, cmap=inferno, alpha=255)\nage_image = set_background(age_image, \"black\")\nage_image"
+ "source": [
+ "# Rasterize the block-group polygons with their PCT_ABOVE_65 attribute.\n",
+ "# Block groups tile the plane (non-overlapping: each pixel is\n # covered by exactly one polygon), so merge='last' is\n",
+ "# equivalent to datashader's ds.mean('PCT_ABOVE_65') here.\n",
+ "age_raster = template.xrs.rasterize(\n",
+ " blockgroup_df, column='PCT_ABOVE_65', merge='last',\n",
+ ")\n",
+ "age_raster.data[~np.isfinite(age_raster.data)] = 0.0\n",
+ "\n",
+ "age_classifed = age_raster.xrs.natural_breaks(num_sample=20000, k=4).where(county_mask)\n",
+ "\n",
+ "age_image = _shade(age_classifed, cmap=inferno, alpha=255)\n",
+ "_show(age_image, bg='black')\n"
+ ]
},
{
"cell_type": "markdown",
@@ -253,7 +354,15 @@
}
},
"outputs": [],
- "source": "pharmacy_deserts = proximity_classifed.xrs.binary([3])\nolder_regions = age_classifed.xrs.binary([3])\ntarget_deserts = (pharmacy_deserts * older_regions).where(county_mask)\ntarget_deserts_img = shade(\n target_deserts, cmap=[\"#333333\", \"fuchsia\"], alpha=255, how=\"linear\"\n)\nset_background(target_deserts_img, \"black\")"
+ "source": [
+ "pharmacy_deserts = proximity_classifed.xrs.binary([3])\n",
+ "older_regions = age_classifed.xrs.binary([3])\n",
+ "target_deserts = (pharmacy_deserts * older_regions).where(county_mask)\n",
+ "target_deserts_img = _shade(\n",
+ " target_deserts, cmap=[\"#333333\", \"fuchsia\"], alpha=255, how=\"linear\"\n",
+ ")\n",
+ "_show(target_deserts_img, bg='black')\n"
+ ]
},
{
"cell_type": "markdown",
@@ -288,11 +397,11 @@
},
"outputs": [],
"source": [
- "from datashader.colors import Set1\n",
- "\n",
- "counties_raster = cvs.polygons(county_df, geometry=\"geometry\", agg=ds.max(\"OBJECTID\"))\n",
- "counties_image = shade(counties_raster, cmap=Set1, alpha=225, how=\"linear\")\n",
- "set_background(counties_image, \"black\")"
+ "counties_raster = template.xrs.rasterize(\n",
+ " county_df, column='OBJECTID', merge='max',\n",
+ ")\n",
+ "counties_image = _shade(counties_raster, cmap=Set1, alpha=225, how=\"linear\")\n",
+ "_show(counties_image, bg='black')\n"
]
},
{
@@ -321,7 +430,19 @@
}
},
"outputs": [],
- "source": "# zones\ncounties_raster.data = counties_raster.data.astype(np.int64)\n\n# values to summarize\ntarget_deserts.data = target_deserts.data.astype(np.int8)\n\n# summary functions\nzonal_funcs = dict(pharmacy_desert_mean=lambda z: z.mean())\n\n# execute summary functions on each zone and take top 10\nresults = target_deserts.xrs.zonal_stats(counties_raster, stats_funcs=zonal_funcs)"
+ "source": [
+ "# zones\n",
+ "counties_raster.data = counties_raster.data.astype(np.int64)\n",
+ "\n",
+ "# values to summarize\n",
+ "target_deserts.data = target_deserts.data.astype(np.int8)\n",
+ "\n",
+ "# summary functions\n",
+ "zonal_funcs = dict(pharmacy_desert_mean=lambda z: z.mean())\n",
+ "\n",
+ "# execute summary functions on each zone and take top 10\n",
+ "results = target_deserts.xrs.zonal_stats(counties_raster, stats_funcs=zonal_funcs)\n"
+ ]
},
{
"cell_type": "markdown",
@@ -345,7 +466,7 @@
"source": [
"cols = [\"pharmacy_desert_mean\", \"NAME\", \"STATE_FIPS\", \"geometry\"]\n",
"final_df = pd.merge(county_df, results, left_on=\"ZONE_ID\", right_index=True)[cols]\n",
- "final_df.nlargest(10, \"pharmacy_desert_mean\")"
+ "final_df.nlargest(10, \"pharmacy_desert_mean\")\n"
]
},
{
@@ -367,7 +488,18 @@
}
},
"outputs": [],
- "source": "counties_raster = cvs.polygons(county_df, geometry=\"geometry\", agg=ds.max(\"OBJECTID\"))\n\ndesert_raster = cvs.polygons(\n final_df.nlargest(75, \"pharmacy_desert_mean\"),\n geometry=\"geometry\",\n agg=ds.mean(\"pharmacy_desert_mean\"),\n)\n\ncounty_mask = cvs.polygons(county_df, geometry=\"geometry\")"
+ "source": [
+ "counties_raster = template.xrs.rasterize(county_df, column='OBJECTID', merge='max')\n",
+ "\n",
+ "# Top-75 counties by pharmacy-desert mean, rasterized with their mean value.\n",
+ "# Counties tile the plane (non-overlapping), so merge='last' matches ds.mean.\n",
+ "desert_raster = template.xrs.rasterize(\n",
+ " final_df.nlargest(75, \"pharmacy_desert_mean\"),\n",
+ " column='pharmacy_desert_mean', merge='last',\n",
+ ")\n",
+ "\n",
+ "county_mask = template.xrs.rasterize(county_df, merge='last')\n"
+ ]
},
{
"cell_type": "code",
@@ -381,7 +513,15 @@
}
},
"outputs": [],
- "source": "img = stack(\n shade(county_mask, cmap=[\"#0f9425\"], alpha=255),\n shade(counties_raster, cmap=[\"green\", \"#ffffff\"], alpha=25),\n shade(desert_raster, cmap=[\"white\", \"#ff0000\"], alpha=255),\n shade(desert_raster.xrs.hillshade(), cmap=[\"green\", \"#ff0000\"], alpha=150),\n)\nset_background(img, \"black\")"
+ "source": [
+ "img = _stack(\n",
+ " _shade(county_mask, cmap=[\"#0f9425\"], alpha=255),\n",
+ " _shade(counties_raster, cmap=[\"green\", \"#ffffff\"], alpha=25),\n",
+ " _shade(desert_raster, cmap=[\"white\", \"#ff0000\"], alpha=255),\n",
+ " _shade(desert_raster.xrs.hillshade(), cmap=[\"green\", \"#ff0000\"], alpha=150),\n",
+ ")\n",
+ "_show(img, bg='black')\n"
+ ]
},
{
"cell_type": "code",
@@ -412,4 +552,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
-}
\ No newline at end of file
+}
diff --git a/examples/xarray-spatial_classification-methods.ipynb b/examples/xarray-spatial_classification-methods.ipynb
index 26cc45192..a2d42244a 100644
--- a/examples/xarray-spatial_classification-methods.ipynb
+++ b/examples/xarray-spatial_classification-methods.ipynb
@@ -47,7 +47,19 @@
},
"outputs": [],
"source": [
- "import xarray as xr\nfrom xrspatial.geotiff import open_geotiff\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = open_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
+ "import matplotlib.pyplot as plt\n",
+ "import xarray as xr\n",
+ "from xrspatial.geotiff import open_geotiff\n",
+ "import xrspatial\n",
+ "\n",
+ "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n",
+ "raster = open_geotiff(file_name, band=0)\n",
+ "raster.name = 'Colorado Elevation Raster'\n",
+ "\n",
+ "xmin, xmax = raster.x.data.min(), raster.x.data.max()\n",
+ "ymin, ymax = raster.y.data.min(), raster.y.data.max()\n",
+ "\n",
+ "xmin, xmax, ymin, ymax\n"
]
},
{
@@ -64,7 +76,7 @@
},
"outputs": [],
"source": [
- "raster.plot.imshow(figsize=(15, 10));"
+ "raster.plot.imshow(figsize=(15, 10));\n"
]
},
{
@@ -104,7 +116,14 @@
}
},
"outputs": [],
- "source": "from datashader.transfer_functions import shade\nimport matplotlib.pyplot as plt\n\nnatural_breaks_agg = raster.xrs.natural_breaks(num_sample=20000, k=15)\n\nshade(natural_breaks_agg, cmap=plt.get_cmap(\"terrain\"), how=\"linear\")"
+ "source": [
+ "natural_breaks_agg = raster.xrs.natural_breaks(num_sample=20000, k=15)\n",
+ "\n",
+ "plt.figure(figsize=(10, 6))\n",
+ "plt.imshow(natural_breaks_agg, cmap=\"terrain\")\n",
+ "plt.colorbar()\n",
+ "plt.show()\n"
+ ]
},
{
"cell_type": "markdown",
@@ -143,7 +162,14 @@
}
},
"outputs": [],
- "source": "equal_interval_agg = raster.xrs.equal_interval(k=15)\n\nshade(equal_interval_agg, cmap=plt.get_cmap(\"terrain\"), how=\"linear\")"
+ "source": [
+ "equal_interval_agg = raster.xrs.equal_interval(k=15)\n",
+ "\n",
+ "plt.figure(figsize=(10, 6))\n",
+ "plt.imshow(equal_interval_agg, cmap=\"terrain\")\n",
+ "plt.colorbar()\n",
+ "plt.show()\n"
+ ]
},
{
"cell_type": "markdown",
@@ -174,7 +200,14 @@
}
},
"outputs": [],
- "source": "quantile_agg = raster.xrs.quantile(k=15)\n\nshade(quantile_agg, cmap=plt.get_cmap(\"terrain\"), how=\"linear\")"
+ "source": [
+ "quantile_agg = raster.xrs.quantile(k=15)\n",
+ "\n",
+ "plt.figure(figsize=(10, 6))\n",
+ "plt.imshow(quantile_agg, cmap=\"terrain\")\n",
+ "plt.colorbar()\n",
+ "plt.show()\n"
+ ]
},
{
"cell_type": "markdown",
@@ -205,7 +238,20 @@
}
},
"outputs": [],
- "source": "bins = range(150, 5200, 350)\nnew_vals = [val if val > 2500 else 0 for val in bins]\n\nreclass_agg = raster.xrs.reclassify(\n bins=bins,\n new_values=new_vals,\n)\n\nshade(reclass_agg, cmap=plt.get_cmap(\"terrain\"), how=\"linear\")"
+ "source": [
+ "bins = range(150, 5200, 350)\n",
+ "new_vals = [val if val > 2500 else 0 for val in bins]\n",
+ "\n",
+ "reclass_agg = raster.xrs.reclassify(\n",
+ " bins=bins,\n",
+ " new_values=new_vals,\n",
+ ")\n",
+ "\n",
+ "plt.figure(figsize=(10, 6))\n",
+ "plt.imshow(reclass_agg, cmap=\"terrain\")\n",
+ "plt.colorbar()\n",
+ "plt.show()\n"
+ ]
},
{
"cell_type": "code",
@@ -237,4 +283,4 @@
},
"nbformat": 4,
"nbformat_minor": 5
-}
\ No newline at end of file
+}
diff --git a/setup.cfg b/setup.cfg
index ef9e1ed2a..7e87a1177 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -52,7 +52,11 @@ doc =
sphinx-panels
sphinx_rtd_theme
examples =
- datashader >= 0.15.0
+ # Rendering (matplotlib) and vector rasterization (geopandas / shapely)
+ # for the example notebooks. datashader is no longer required.
+ matplotlib
+ geopandas
+ shapely
plot =
# Optional for the .xrs.plot accessor helpers. Every matplotlib
# import in the package is lazy, so the compute functions work
diff --git a/xrspatial/terrain.py b/xrspatial/terrain.py
index b6f4a0ce2..17ac10509 100644
--- a/xrspatial/terrain.py
+++ b/xrspatial/terrain.py
@@ -885,7 +885,7 @@ def generate_terrain(agg: xr.DataArray,
if 'x' in agg.coords and 'y' in agg.coords:
coords = {'y': agg.coords['y'], 'x': agg.coords['x']}
else:
- # pixel-center coordinates (matches datashader Canvas convention)
+ # pixel-center coordinates (standard pixel-center grid convention)
dx = (x_range[1] - x_range[0]) / width
dy = (y_range[1] - y_range[0]) / height
xs = np.linspace(x_range[0] + dx / 2, x_range[1] - dx / 2, width)
diff --git a/xrspatial/tests/test_utils.py b/xrspatial/tests/test_utils.py
index f8ac187ea..84d415a94 100644
--- a/xrspatial/tests/test_utils.py
+++ b/xrspatial/tests/test_utils.py
@@ -301,3 +301,25 @@ def test_slope_irregular_coords_warns():
)
with pytest.warns(UserWarning, match="'x' coordinate is not evenly spaced"):
slope(raster)
+
+
+def test_canvas_like_removed_3451_3640():
+ """canvas_like() was the only library-source datashader import.
+
+ It was removed in #3451 / #3640 (decouple datashader). Guard against
+ accidental re-introduction: the attribute must be absent and
+ `xrspatial.utils` must not import datashader at module load time.
+ """
+ import xrspatial.utils as u
+
+ assert not hasattr(u, "canvas_like"), (
+ "canvas_like was removed in #3451; do not re-add it (it pulled "
+ "datashader into the library source)."
+ )
+
+ # No datashader import should be referenced from the utils module.
+ import sys
+
+ assert "datashader" not in sys.modules, (
+ "importing xrspatial.utils should not load datashader."
+ )
diff --git a/xrspatial/utils.py b/xrspatial/utils.py
index 44d7e5a5f..041362b68 100644
--- a/xrspatial/utils.py
+++ b/xrspatial/utils.py
@@ -791,84 +791,6 @@ def bands_to_img(r, g, b, nodata=1):
return Image.fromarray(data, "RGBA")
-def canvas_like(
- raster,
- width=512,
- height=None,
- x_range=None,
- y_range=None,
- **kwargs
-):
-
- """
- Resample a xarray.DataArray by canvas width and bounds.
- Height of the resampled raster is implied from the canvas width
- using aspect ratio of original raster.
-
- This function uses of datashader.Canvas.raster internally.
- Most of the docstrings are copied from Datashader.
-
- Handles 2D or 3D xarray.DataArray, assuming that the last two
- array dimensions are the y-axis and x-axis that are to be
- resampled. If a 3D array is supplied a layer may be specified
- to resample to select the layer along the first dimension to
- resample.
-
- Parameters
- ----------
- raster : xarray.DataArray
- 2D or 3D labeled data array.
- layer : float, optional
- For a 3D array, value along the z dimension.
- width : int, default=512
- Width of the output aggregate in pixels.
- height : int, default=None
- Height of the output aggregate in pixels.
- If not provided, height will be implied from `width`
- using aspect ratio of input raster.
- x_range : tuple of int, optional
- A tuple representing the bounds inclusive space ``[min, max]``
- along the x-axis.
- y_range : tuple of int, optional
- A tuple representing the bounds inclusive space ``[min, max]``
- along the y-axis.
-
- References
- ----------
- - https://datashader.org/_modules/datashader/core.html#Canvas
- """
-
- # get ranges
- if x_range is None:
- x_range = (
- raster.coords["x"].min().item(),
- raster.coords["x"].max().item()
- )
- if y_range is None:
- y_range = (
- raster.coords["y"].min().item(),
- raster.coords["y"].max().item()
- )
-
- if height is None:
- # set width and height
- height = height_implied_by_aspect_ratio(width, x_range, y_range)
-
- try:
- import datashader as ds
- except ImportError:
- raise ImportError(
- "canvas_like requires datashader: pip install datashader"
- )
-
- cvs = ds.Canvas(
- plot_width=width, plot_height=height, x_range=x_range, y_range=y_range
- )
- out = cvs.raster(raster, **kwargs)
-
- return out
-
-
def _hex_to_rgb(c):
"""Convert a hex color string (e.g. '#ff0000' or 'ff0000') to (r, g, b)."""
c = c.lstrip('#')
diff --git a/xrspatial/zonal.py b/xrspatial/zonal.py
index 7167de67e..6e92f96be 100644
--- a/xrspatial/zonal.py
+++ b/xrspatial/zonal.py
@@ -2112,10 +2112,9 @@ def suggest_zonal_canvas(
--------
.. sourcecode:: python
- >>> # Imports (datashader is optional: pip install datashader)
- >>> from spatialpandas import GeoDataFrame
+ >>> import numpy as np
+ >>> import xarray as xr
>>> import geopandas as gpd
- >>> import datashader as ds
>>> from xrspatial.zonal import suggest_zonal_canvas
>>> df = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
@@ -2141,10 +2140,14 @@ def suggest_zonal_canvas(
)
>>> height, width
(1537, 2376)
- >>> cvs = ds.Canvas(x_range=x_range, y_range=y_range,
- >>> plot_height=height, plot_width=width)
- >>> spatial_df = GeoDataFrame(df, geometry='geometry')
- >>> agg = cvs.polygons(spatial_df, 'geometry', agg=ds.max('id'))
+ >>> # Build a template grid and rasterize the polygons onto it with
+ >>> # the .xrs.rasterize accessor (no datashader needed).
+ >>> xs = np.linspace(x_range[0], x_range[1], width)
+ >>> ys = np.linspace(y_range[0], y_range[1], height)
+ >>> template = xr.DataArray(
+ ... np.full((height, width), np.nan),
+ ... coords={'y': ys, 'x': xs}, dims=['y', 'x'])
+ >>> agg = template.xrs.rasterize(df, column='id', merge='max')
>>> min_poly_id = df.area.argmin()
>>> actual_min_pixels = len(np.where(agg.data==min_poly_id)[0])
>>> actual_min_pixels
| \n", - " | zone | \n", - "elevation_change | \n", - "elevation_min | \n", - "elevation_max | \n", - "
|---|---|---|---|---|
| 0 | \n", - "0 | \n", - "4000.000000 | \n", - "0.000000 | \n", - "4000.000000 | \n", - "
| 1 | \n", - "11 | \n", - "16.527466 | \n", - "1435.497249 | \n", - "1452.024715 | \n", - "
| 2 | \n", - "70 | \n", - "0.000000 | \n", - "1435.677678 | \n", - "1435.677678 | \n", - "