Surface#
With the Surface tools, you can quantify and visualize a terrain landform represented by a digital elevation model.
Starting with a raster elevation surface that represented as an Xarray DataArray (or an Xarray Dataset containing multiple elevation variables), these tools help you in identifying some specific patterns that were not readily apparent in the original surface. When a DataArray is passed, the return is a DataArray. When a Dataset is passed, the function is applied to each variable independently and the return is a Dataset.
Hillshade: Creates a shaded relief from a surface raster by considering the illumination source angle and shadows.
Slope: Identifies the slope from each cell of a raster.
Curvature: Calculates the curvature of a raster surface.
Aspect: Derives the aspect from each cell of a raster surface.
Viewshed: Determines visible locations in the input raster surface from a viewpoint with some optional observer features.
Importing Packages#
[1]:
import numpy as np
import xarray as xr
import xrspatial
[2]:
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
import numpy as np
def _shade(arr, cmap=None, alpha=None, min_alpha=None, span=None, how='linear'):
arr = np.asarray(arr, dtype=np.float64)
finite = np.isfinite(arr)
if span is not None:
lo, hi = span
elif finite.any():
lo = np.nanmin(arr); hi = np.nanmax(arr)
else:
lo = hi = 0.0
norm = np.zeros_like(arr, dtype=np.float32)
if hi > lo:
norm = np.where(finite, (arr - lo) / (hi - lo), 0.0).astype(np.float32)
if isinstance(cmap, (list, tuple)):
cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
elif cmap is None:
cmap = plt.get_cmap('terrain')
rgba = cmap(norm, bytes=True)
a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
if min_alpha is not None:
a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
rgba[..., 3] = a
return rgba
def _stack(*layers):
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], base[..., 3])
base = np.concatenate([rgb, out_a[..., None]], axis=-1)
return np.clip(base, 0, 255).astype(np.uint8)
def _show(rgba, bg=None, figsize=(8, 6)):
plt.figure(figsize=figsize)
plt.imshow(rgba)
plt.axis('off')
plt.show()
Elevation = LinearSegmentedColormap.from_list(
'elevation',
[(0.0, (0.0, 0.3, 0.6)), (0.2, (0.2, 0.6, 0.85)), (0.3, (0.55, 0.78, 0.92)),
(0.32, (0.4, 0.73, 0.45)), (0.45, (0.3, 0.55, 0.25)), (0.6, (0.5, 0.5, 0.3)),
(0.75, (0.65, 0.55, 0.4)), (0.9, (0.8, 0.8, 0.78)), (1.0, (1.0, 1.0, 1.0))],
N=256)
Generate Terrain Data#
The rest of the geo-related functions focus on raster data, i.e. data that’s been aggregated into the row-column grid of cells in a raster image. Datashader’s Canvas object provides a convenient frame to set up a new raster, so we’ll use that with our generate_terrain function to generate some fake terrain as an elevation raster. Once we have that, we’ll use matplotlib’s shade for easy visualization.
[3]:
from xrspatial import generate_terrain
W = 800
H = 600
terrain = xr.DataArray(np.zeros((H, W)))
terrain = generate_terrain(terrain)
_show(_shade(terrain, cmap=["black", "white"], how="linear"))
The grayscale values in the image above show elevation, scaled linearly in black-to-white color intensity (with the large black areas indicating low elevation). This shows the data, but it would look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops. Let’s try the Elevation colormap we imported above:
[4]:
_show(_shade(terrain, cmap=Elevation, how="linear"))
Hillshade#
Hillshade is a technique used to visualize terrain as shaded relief by illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which can be calculated from slope and aspect.
Let’s apply Hillshade to our terrain and visualize the result with shade.
[5]:
from xrspatial import hillshade
illuminated = hillshade(terrain)
hillshade_gray_white = _shade(
illuminated, cmap=["gray", "white"], alpha=255, how="linear"
)
hillshade_gray_white
[5]:
array([[[128, 128, 128, 0],
[128, 128, 128, 0],
[128, 128, 128, 0],
...,
[128, 128, 128, 0],
[128, 128, 128, 0],
[128, 128, 128, 0]],
[[128, 128, 128, 0],
[187, 187, 187, 255],
[187, 187, 187, 255],
...,
[187, 187, 187, 255],
[187, 187, 187, 255],
[128, 128, 128, 0]],
[[128, 128, 128, 0],
[187, 187, 187, 255],
[187, 187, 187, 255],
...,
[187, 187, 187, 255],
[187, 187, 187, 255],
[128, 128, 128, 0]],
...,
[[128, 128, 128, 0],
[187, 187, 187, 255],
[187, 187, 187, 255],
...,
[187, 187, 187, 255],
[187, 187, 187, 255],
[128, 128, 128, 0]],
[[128, 128, 128, 0],
[187, 187, 187, 255],
[187, 187, 187, 255],
...,
[187, 187, 187, 255],
[187, 187, 187, 255],
[128, 128, 128, 0]],
[[128, 128, 128, 0],
[128, 128, 128, 0],
[128, 128, 128, 0],
...,
[128, 128, 128, 0],
[128, 128, 128, 0],
[128, 128, 128, 0]]], shape=(600, 800, 4), dtype=uint8)
Applying hillshade reveals a lot of detail in the 3D shape of the terrain.
To add even more detail, we can add the Elevation colormapped terrain from earlier and combine it with the hillshade terrain using matplotlib’s stack function.
[6]:
terrain_elevation = _shade(terrain, cmap=Elevation, alpha=128, how="linear")
_show(_stack(hillshade_gray_white, terrain_elevation))
Slope#
Slope is the inclination of a surface. In geography, slope is the amount of change in elevation for an area in a terrain relative to its surroundings.
Xarray-spatial’s slope function returns the slope at each cell in degrees. Because Xarray-spatial is integrated with Xarray and Numpy, we can apply standard Numpy filters. For example, we can highlight only slopes in the avalanche risk range of 25 - 50 degrees. (Note the use of risky.data since these are DataArrays). Stacking the resulting raster with the hillshaded and plain terrain ones from above gives an image with areas of avalanche risk neatly highlighted.
[7]:
from xrspatial import slope
risky = slope(terrain)
risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)
_show(_stack(
_shade(terrain, cmap=["black", "white"], how="linear"),
_shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
_shade(risky, cmap="red", how="linear", alpha=200),
))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[7], line 9
5
6 _show(_stack(
7 _shade(terrain, cmap=["black", "white"], how="linear"),
8 _shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
----> 9 _shade(risky, cmap="red", how="linear", alpha=200),
10 ))
Cell In[2], line 21, in _shade(arr, cmap, alpha, min_alpha, span, how)
17 if isinstance(cmap, (list, tuple)):
18 cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
19 elif cmap is None:
20 cmap = plt.get_cmap('terrain')
---> 21 rgba = cmap(norm, bytes=True)
22 a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
23 if min_alpha is not None:
24 a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
TypeError: 'str' object is not callable
Curvature#
Curvature is the second derivative of a surface’s elevation, or the slope-of-the-slope; in other words, how fast the slope is increasing or decreasing as we move along a surface.
A positive curvature means the surface is curving up (upwardly convex) at that cell.
A negative curvature means the surface is curving down (downwardly convex) at that cell.
A curvature of 0 means the surface is striaght and constant in whatever angle it’s sloped towards.
Let’s generate a terrain with an appropriate z-factor and apply the curvature function to it. Then, we can apply some Numpy filtering (remember, we have access to all those functions) to highlight steeper and gentler curves in the slopes. Stacking these with the hillshaded and plain terrains gives us a fuller picture of the slopes.
[8]:
from xrspatial import curvature
terrain_z_one = xr.DataArray(np.zeros((H, W)))
terrain_z_one = generate_terrain(terrain_z_one, zfactor=1)
curv = curvature(terrain_z_one)
curv_hi, curv_low = curv.copy(), curv.copy()
curv_hi.data = np.where(np.logical_and(curv_hi.data > 1, curv_hi.data < 4), 1, np.nan)
curv_low.data = np.where(
np.logical_and(curv_low.data > 0.5, curv_low.data < 1), 1, np.nan
)
_show(_stack(
_shade(terrain, cmap=["black", "white"], how="linear"),
_shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
_shade(curv_hi, cmap="red", how="log", alpha=200),
_shade(curv_low, cmap="green", how="log", alpha=200),
))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[8], line 15
11
12 _show(_stack(
13 _shade(terrain, cmap=["black", "white"], how="linear"),
14 _shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
---> 15 _shade(curv_hi, cmap="red", how="log", alpha=200),
16 _shade(curv_low, cmap="green", how="log", alpha=200),
17 ))
Cell In[2], line 21, in _shade(arr, cmap, alpha, min_alpha, span, how)
17 if isinstance(cmap, (list, tuple)):
18 cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
19 elif cmap is None:
20 cmap = plt.get_cmap('terrain')
---> 21 rgba = cmap(norm, bytes=True)
22 a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
23 if min_alpha is not None:
24 a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
TypeError: 'str' object is not callable
Aspect#
Aspect is the orientation of a slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.
The Xarray-spatial aspect function returns the aspect in degrees for each cell in an elevation terrain.
We can apply aspect to our terrain, then use Numpy to filter out only slopes facing close to North. Then, we can stack that with the hillshaded and plain terrains. (Note: the printout images are from a North point-of-view.)
[9]:
from xrspatial import aspect
north_faces = aspect(terrain)
north_faces.data = np.where(
np.logical_or(north_faces.data > 350, north_faces.data < 10), 1, np.nan
)
_show(_stack(
_shade(terrain, cmap=["black", "white"], how="linear"),
_shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
_shade(north_faces, cmap=["aqua"], how="linear", alpha=100),
))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[9], line 11
7
8 _show(_stack(
9 _shade(terrain, cmap=["black", "white"], how="linear"),
10 _shade(illuminated, cmap=["black", "white"], how="linear", alpha=128),
---> 11 _shade(north_faces, cmap=["aqua"], how="linear", alpha=100),
12 ))
Cell In[2], line 21, in _shade(arr, cmap, alpha, min_alpha, span, how)
17 if isinstance(cmap, (list, tuple)):
18 cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
19 elif cmap is None:
20 cmap = plt.get_cmap('terrain')
---> 21 rgba = cmap(norm, bytes=True)
22 a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
23 if min_alpha is not None:
24 a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/colors.py:779, in Colormap.__call__(self, X, alpha, bytes)
756 def __call__(self, X, alpha=None, bytes=False):
757 r"""
758 Parameters
759 ----------
(...) 777 RGBA values with a shape of ``X.shape + (4, )``.
778 """
--> 779 rgba, mask = self._get_rgba_and_mask(X, alpha=alpha, bytes=bytes)
780 if not np.iterable(X):
781 rgba = tuple(rgba)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/colors.py:809, in Colormap._get_rgba_and_mask(self, X, alpha, bytes)
784 def _get_rgba_and_mask(self, X, alpha=None, bytes=False):
785 r"""
786 Parameters
787 ----------
(...) 807 Boolean array with True where the input is ``np.nan`` or masked.
808 """
--> 809 self._ensure_inited()
811 xa = np.array(X, copy=True)
812 if not xa.dtype.isnative:
813 # Native byteorder is faster.
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/colors.py:982, in Colormap._ensure_inited(self)
980 def _ensure_inited(self):
981 if not self._isinit:
--> 982 self._init()
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/colors.py:1167, in LinearSegmentedColormap._init(self)
1164 def _init(self):
1165 # Assemble the LUT first in a local variable in case of parallel threads
1166 lut = np.ones((self.N + 3, 4), float)
-> 1167 lut[:-3, 0] = _create_lookup_table(
1168 self.N, self._segmentdata['red'], self._gamma)
1169 lut[:-3, 1] = _create_lookup_table(
1170 self.N, self._segmentdata['green'], self._gamma)
1171 lut[:-3, 2] = _create_lookup_table(
1172 self.N, self._segmentdata['blue'], self._gamma)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/colors.py:682, in _create_lookup_table(N, data, gamma)
679 y1 = adata[:, 2]
681 if x[0] != 0. or x[-1] != 1.0:
--> 682 raise ValueError(
683 "data mapping points must start with x=0 and end with x=1")
684 if (np.diff(x) < 0).any():
685 raise ValueError("data mapping points must have x in increasing order")
ValueError: data mapping points must start with x=0 and end with x=1
Viewshed#
The xrspatial.viewshed function operates on a given aggregate to calculate the viewshed (the visible cells in the raster) for a given viewpoint, or observer location.
The visibility model is as follows: Two cells are visible to each other if the line of sight that connects their centers is not blocked at any point by another part of the terrain. If the line of sight does not pass through the cell center, elevation is determined using bilinear interpolation.
Simple Viewshed Example#
The example below creates a raster aggregate from a 2d normal distribution.
To calculate the viewshed, we need an observer location so we’ll set up an aggregate for that as well.
Then, we can visualize all of that with hillshade, shade, and stack.
The observer location is indicated by the orange point in the upper-left of the plot.
[10]:
from xrspatial import viewshed
import pandas as pd
import xarray as xr
OBSERVER_X = -12.5
OBSERVER_Y = 10
W = 400
H = 300
x_range = (-20, 20)
y_range = (-20, 20)
# Rasterize a 2D normal distribution into a density grid using a plain
# numpy 2D histogram (replaces ds.Canvas().points(...) count aggregation).
normal_df = pd.DataFrame(
{"x": np.random.normal(0.5, 1, 10_000_000),
"y": np.random.normal(0.5, 1, 10_000_000)}
)
xs = np.linspace(x_range[0], x_range[1], W)
ys = np.linspace(y_range[1], y_range[0], H)
counts, _, _ = np.histogram2d(normal_df['x'], normal_df['y'],
bins=[W, H],
range=[[x_range[0], x_range[1]],
[y_range[0], y_range[1]]])
# histogram2d returns shape (W, H) indexed by x then y; orient to (H, W)
normal_agg = xr.DataArray(counts.T.astype('float64'),
coords={'y': ys, 'x': xs}, dims=['y', 'x'])
normal_shaded = _shade(normal_agg)
# Single observer point via the .xrs.rasterize accessor.
import geopandas as gpd
from shapely.geometry import Point
template = xr.DataArray(np.full((H, W), np.nan),
coords={'y': ys, 'x': xs}, dims=['y', 'x'])
observer_gdf = gpd.GeoDataFrame(
{'id': [1]}, geometry=[Point(OBSERVER_X, OBSERVER_Y)], crs='EPSG:4326')
observer_agg = template.xrs.rasterize(observer_gdf, column='id', merge='last')
observer_shaded = _shade(observer_agg, cmap=['orange', 'orange'], min_alpha=255)
normal_illuminated = hillshade(normal_agg)
normal_illuminated_shaded = _shade(
normal_illuminated, cmap=["black", "white"], alpha=128, how="linear"
)
_show(_stack(normal_illuminated_shaded, observer_shaded))
Calculate viewshed using the observer location#
Now we can apply viewshed to the normal_agg, with the observer_agg for the viewpoint. We can then visualize it and stack it with the hillshade and observer rasters.
[11]:
# Will take some time to run...
%time view = viewshed(normal_agg, x=OBSERVER_X, y=OBSERVER_Y)
view_shaded = _shade(view, cmap=["white", "red"], alpha=128, how="linear")
_show(_stack(normal_illuminated_shaded, observer_shaded, view_shaded))
CPU times: user 14.9 s, sys: 1.22 s, total: 16.2 s
Wall time: 16.2 s
As you can see, the image highlights in red all points visible from the observer location marked with the orange dot. As one might expect, the areas behind the normal distribution mountain are blocked from the viewer.
Viewshed on Terrain#
Now we can try using viewshed on our more complicated terrain.
We’ll set up our terrain aggregate and apply hillshade and shade for easy visualization.
We’ll also set up an observer location aggregate, setting the location to the center, at (x, y) = (0, 0).
[12]:
from xrspatial import viewshed
import geopandas as gpd
from shapely.geometry import Point
x_range = (-20e6, 20e6)
y_range = (-20e6, 20e6)
terrain = xr.DataArray(np.zeros((H, W)))
terrain = generate_terrain(terrain, x_range=x_range, y_range=y_range)
terrain_shaded = _shade(terrain, cmap=Elevation, alpha=128, how="linear")
illuminated = hillshade(terrain)
OBSERVER_X = 0.0
OBSERVER_Y = 0.0
# Observer point via the .xrs.rasterize accessor (replaces ds.Canvas().points).
xs = np.linspace(x_range[0], x_range[1], W)
ys = np.linspace(y_range[1], y_range[0], H)
template = xr.DataArray(np.full((H, W), np.nan),
coords={'y': ys, 'x': xs}, dims=['y', 'x'])
observer_gdf = gpd.GeoDataFrame(
{'id': [1]}, geometry=[Point(OBSERVER_X, OBSERVER_Y)], crs='EPSG:4326')
observer_agg = template.xrs.rasterize(observer_gdf, column='id', merge='last')
observer_shaded = _shade(observer_agg, cmap=['orange', 'orange'], min_alpha=255)
_show(_stack(
_show(_shade(illuminated, cmap=["black", "white"], alpha=128, how="linear"),
terrain_shaded,
observer_shaded,
)))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[12], line 27
23 observer_agg = template.xrs.rasterize(observer_gdf, column='id', merge='last')
24 observer_shaded = _shade(observer_agg, cmap=['orange', 'orange'], min_alpha=255)
25
26 _show(_stack(
---> 27 _show(_shade(illuminated, cmap=["black", "white"], alpha=128, how="linear"),
28 terrain_shaded,
29 observer_shaded,
30 )))
Cell In[2], line 38, in _show(rgba, bg, figsize)
37 def _show(rgba, bg=None, figsize=(8, 6)):
---> 38 plt.figure(figsize=figsize)
39 plt.imshow(rgba)
40 plt.axis('off')
41 plt.show()
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/pyplot.py:1096, in figure(num, figsize, dpi, facecolor, edgecolor, frameon, FigureClass, clear, **kwargs)
1086 if len(allnums) == max_open_warning >= 1:
1087 _api.warn_external(
1088 f"More than {max_open_warning} figures have been opened. "
1089 f"Figures created through the pyplot interface "
(...) 1093 f"Consider using `matplotlib.pyplot.close()`.",
1094 RuntimeWarning)
-> 1096 manager = new_figure_manager(
1097 num, figsize=figsize, dpi=dpi,
1098 facecolor=facecolor, edgecolor=edgecolor, frameon=frameon,
1099 FigureClass=FigureClass, **kwargs)
1100 fig = manager.canvas.figure
1101 if fig_label:
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/pyplot.py:576, in new_figure_manager(*args, **kwargs)
574 """Create a new figure manager instance."""
575 _warn_if_gui_out_of_main_thread()
--> 576 return _get_backend_mod().new_figure_manager(*args, **kwargs)
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib_inline/backend_inline.py:26, in new_figure_manager(num, FigureClass, *args, **kwargs)
20 def new_figure_manager(num, *args, FigureClass=Figure, **kwargs):
21 """
22 Return a new figure manager for a new figure instance.
23
24 This function is part of the API expected by Matplotlib backends.
25 """
---> 26 return new_figure_manager_given_figure(num, FigureClass(*args, **kwargs))
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/figure.py:2632, in Figure.__init__(self, figsize, dpi, facecolor, edgecolor, linewidth, frameon, subplotpars, tight_layout, constrained_layout, layout, **kwargs)
2629 edgecolor = mpl._val_or_rc(edgecolor, 'figure.edgecolor')
2630 frameon = mpl._val_or_rc(frameon, 'figure.frameon')
-> 2632 figsize = _parse_figsize(figsize, dpi)
2634 if not np.isfinite(figsize).all() or (np.array(figsize) < 0).any():
2635 raise ValueError('figure size must be positive finite not '
2636 f'{figsize}')
File ~/checkouts/readthedocs.org/user_builds/xarray-spatial/envs/stable/lib/python3.12/site-packages/matplotlib/figure.py:3794, in _parse_figsize(figsize, dpi)
3789 raise ValueError(
3790 f"Invalid unit {unit!r} in 'figsize'; "
3791 "supported units are 'in', 'cm', 'px'"
3792 )
3793 else:
-> 3794 raise ValueError(
3795 "Invalid figsize format, expected (x, y) or (x, y, unit) but got "
3796 f"{figsize!r}"
3797 )
3799 if x is None and y is None:
3800 raise ValueError(
3801 "figsize=(None, None) is invalid; at least one of width or "
3802 "height must be provided")
ValueError: Invalid figsize format, expected (x, y) or (x, y, unit) but got array([[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]],
[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]],
[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]],
...,
[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]],
[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]],
[[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0],
...,
[255, 165, 0, 0],
[255, 165, 0, 0],
[255, 165, 0, 0]]], shape=(300, 400, 4), dtype=uint8)
Now we can apply viewshed.
Notice the use of the
observer_elevargument, which is the height of the observer above the terrain.
[13]:
%time view = viewshed(terrain, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=100)
view_shaded = _shade(view, cmap="fuchsia", how="linear")
_show(_stack(
_shade(illuminated, cmap=["black", "white"], alpha=128, how="linear"),
terrain_shaded,
view_shaded,
observer_shaded,
))
CPU times: user 2.12 s, sys: 4.93 ms, total: 2.12 s
Wall time: 2.12 s
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[13], line 3
1 get_ipython().run_line_magic('time', 'view = viewshed(terrain, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=100)')
2
----> 3 view_shaded = _shade(view, cmap="fuchsia", how="linear")
4 _show(_stack(
5 _shade(illuminated, cmap=["black", "white"], alpha=128, how="linear"),
6 terrain_shaded,
Cell In[2], line 21, in _shade(arr, cmap, alpha, min_alpha, span, how)
17 if isinstance(cmap, (list, tuple)):
18 cmap = LinearSegmentedColormap.from_list('c', list(cmap), N=256)
19 elif cmap is None:
20 cmap = plt.get_cmap('terrain')
---> 21 rgba = cmap(norm, bytes=True)
22 a = np.where(finite, 255 if alpha is None else alpha, 0).astype(np.uint8)
23 if min_alpha is not None:
24 a = np.where(finite & (a < min_alpha), min_alpha, a).astype(np.uint8)
TypeError: 'str' object is not callable
The fuchsia areas are those visible to an observer of the given height at the indicated orange location.
References#
An overview of the Surface toolset: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/an-overview-of-the-surface-tools.htm
Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), p. 406.
Making Maps with Noise Functions: https://www.redblobgames.com/maps/terrain-from-noise/
How Aspect Works: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-aspect-works.htm#ESRI_SECTION1_4198691F8852475A9F4BC71246579FAA