Skip to content

Commit 58475ad

Browse files
brendancolbrendancol
authored
added minimap to explore() mode; added voxelate to mesh_utils (#48)
* added minimap to explore() mode; added voxelate to mesh_utils * changed mesh type names * exposed more of the interactive viewer controls and fixed x-axis flip bug * made trinidad interactive viewer larger * added fix for blockdim on voxelate * way too big a commit in way too big a PR * clamping the observer to terrain...no below terrain for now... * added in some new examples cities * incremented version
1 parent ada1318 commit 58475ad

28 files changed

Lines changed: 9602 additions & 1234 deletions

‎.gitignore‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,3 +37,7 @@ build/
3737
dist/
3838
rtxpy.egg-info/
3939
examples/.ipynb_checkpoints/
40+
*.npz
41+
*.png
42+
*.tif
43+
*.geojson

‎README.md‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,28 @@ Ray tracing using CUDA, accessible from Python.
66

77
*Real-time viewshed analysis with GPU-accelerated ray tracing. Green areas are visible from the observer position (blue dot). Run `python examples/playground.py` to try it interactively.*
88

9+
## Quick Start
10+
11+
```python
12+
import rtxpy # registers the .rtx xarray accessor
13+
import rioxarray
14+
15+
# Load a GeoTIFF DEM as an xarray DataArray
16+
dem = rioxarray.open_rasterio('elevation.tif').squeeze()
17+
18+
# Move data to GPU
19+
dem = dem.rtx.to_cupy()
20+
21+
# Compute hillshade with ray-traced shadows
22+
hillshade = dem.rtx.hillshade(shadows=True)
23+
24+
# Compute viewshed from an observer location (pixel coordinates)
25+
viewshed = dem.rtx.viewshed(x=500, y=300, observer_elev=2)
26+
27+
# Launch interactive 3D terrain explorer
28+
dem.rtx.explore()
29+
```
30+
931
## Prerequisites
1032

1133
- NVIDIA GPU with RTX support (Maxwell architecture or newer)

‎examples/capetown.py‎

Lines changed: 292 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,292 @@
1+
"""Interactive playground for Cape Town, South Africa.
2+
3+
Explore the terrain of Cape Town using GPU-accelerated ray tracing.
4+
Elevation data is sourced from the Copernicus GLO-30 DEM (30 m).
5+
6+
Builds an xr.Dataset with elevation, slope, aspect, and quantile layers.
7+
Press G to cycle between layers. Satellite tiles are draped on the terrain
8+
automatically — press U to toggle tile overlay on/off.
9+
10+
Requirements:
11+
pip install rtxpy[all] matplotlib xarray rioxarray requests pyproj Pillow
12+
"""
13+
14+
import numpy as np
15+
import xarray as xr
16+
17+
from xrspatial import slope, aspect, quantile
18+
from pathlib import Path
19+
20+
import warnings
21+
22+
from rtxpy import fetch_dem, fetch_buildings, fetch_roads, fetch_water, fetch_firms
23+
import rtxpy
24+
25+
# Water feature classification
26+
_MAJOR_WATER = {'river', 'canal'}
27+
_MINOR_WATER = {'stream', 'drain', 'ditch'}
28+
29+
# Cape Town bounding box (lon_min, lat_min, lon_max, lat_max)
30+
BOUNDS = (18.3, -34.2, 18.7, -33.8)
31+
32+
33+
def load_terrain():
34+
"""Load Cape Town terrain data, downloading if necessary."""
35+
dem_path = Path(__file__).parent / "capetown_dem.tif"
36+
37+
terrain = fetch_dem(
38+
bounds=BOUNDS,
39+
output_path=dem_path,
40+
source='copernicus',
41+
crs='EPSG:32734', # UTM zone 34S
42+
)
43+
44+
# Scale down elevation for visualization (optional)
45+
terrain.data = terrain.data * 0.025
46+
47+
# Ensure contiguous array before GPU transfer
48+
terrain.data = np.ascontiguousarray(terrain.data)
49+
50+
# Get stats before GPU transfer (nanmin/nanmax to skip NaN ocean pixels)
51+
elev_min = float(np.nanmin(terrain.data))
52+
elev_max = float(np.nanmax(terrain.data))
53+
54+
# Convert to cupy for GPU processing using the accessor
55+
terrain = terrain.rtx.to_cupy()
56+
57+
print(f"Terrain loaded: {terrain.shape}, elevation range: "
58+
f"{elev_min:.0f}m to {elev_max:.0f}m (scaled)")
59+
60+
return terrain
61+
62+
63+
if __name__ == "__main__":
64+
# Load terrain data (downloads if needed)
65+
terrain = load_terrain()
66+
67+
print("\nControls:")
68+
print(" W/S/A/D or Arrow keys: Move camera")
69+
print(" Q/E or Page Up/Down: Move up/down")
70+
print(" I/J/K/L: Look around")
71+
print(" +/-: Adjust movement speed")
72+
print(" G: Cycle overlay layers")
73+
print(" O: Place observer (for viewshed)")
74+
print(" V: Toggle viewshed (teal glow)")
75+
print(" [/]: Adjust observer height")
76+
print(" T: Toggle shadows")
77+
print(" C: Cycle colormap")
78+
print(" U: Toggle tile overlay")
79+
print(" F: Screenshot")
80+
print(" H: Toggle help overlay")
81+
print(" X: Exit\n")
82+
83+
# Build Dataset with derived layers
84+
print("Building Dataset with terrain analysis layers...")
85+
ds = xr.Dataset({
86+
'elevation': terrain.rename(None),
87+
'slope': slope(terrain),
88+
'aspect': aspect(terrain),
89+
'quantile': quantile(terrain),
90+
})
91+
print(ds)
92+
93+
# Drape satellite tiles on terrain (reprojected to match DEM CRS)
94+
print("Loading satellite tiles...")
95+
ds.rtx.place_tiles('satellite', z='elevation')
96+
97+
# --- Microsoft Global Building Footprints --------------------------------
98+
try:
99+
bldg_cache = Path(__file__).parent / "capetown_buildings.geojson"
100+
bldg_data = fetch_buildings(
101+
bounds=BOUNDS,
102+
cache_path=bldg_cache,
103+
)
104+
105+
# Scale building heights to match the 0.025× terrain elevation.
106+
elev_scale = 0.025
107+
default_height_m = 8.0
108+
for feat in bldg_data.get("features", []):
109+
props = feat.get("properties", {})
110+
h = props.get("height", -1)
111+
if not isinstance(h, (int, float)) or h <= 0:
112+
h = default_height_m
113+
props["height"] = h * elev_scale
114+
115+
mesh_cache_path = Path(__file__).parent / "capetown_buildings_mesh.npz"
116+
with warnings.catch_warnings():
117+
warnings.filterwarnings("ignore", message="place_geojson called before")
118+
bldg_info = ds.rtx.place_geojson(
119+
bldg_data,
120+
z='elevation',
121+
height=default_height_m * elev_scale,
122+
height_field='height',
123+
geometry_id='building',
124+
densify=False,
125+
merge=True,
126+
extrude=True,
127+
mesh_cache=mesh_cache_path,
128+
)
129+
print(f"Placed {bldg_info['geometries']} building footprint geometries")
130+
except ImportError as e:
131+
print(f"Skipping buildings: {e}")
132+
133+
# --- OpenStreetMap roads ------------------------------------------------
134+
try:
135+
# Major roads: motorways, trunk, primary, secondary
136+
major_cache = Path(__file__).parent / "capetown_roads_major.geojson"
137+
major_roads = fetch_roads(
138+
bounds=BOUNDS,
139+
road_type='major',
140+
cache_path=major_cache,
141+
)
142+
if major_roads.get('features'):
143+
with warnings.catch_warnings():
144+
warnings.filterwarnings("ignore", message="place_geojson called before")
145+
info = ds.rtx.place_geojson(
146+
major_roads, z='elevation', height=1,
147+
label_field='name', geometry_id='road_major',
148+
color=(0.10, 0.10, 0.10),
149+
densify=False,
150+
merge=True,
151+
mesh_cache=Path(__file__).parent / "capetown_roads_major_mesh.npz",
152+
)
153+
print(f"Placed {info['geometries']} major road geometries")
154+
155+
# Minor roads: tertiary, residential, service
156+
minor_cache = Path(__file__).parent / "capetown_roads_minor.geojson"
157+
minor_roads = fetch_roads(
158+
bounds=BOUNDS,
159+
road_type='minor',
160+
cache_path=minor_cache,
161+
)
162+
if minor_roads.get('features'):
163+
with warnings.catch_warnings():
164+
warnings.filterwarnings("ignore", message="place_geojson called before")
165+
info = ds.rtx.place_geojson(
166+
minor_roads, z='elevation', height=1,
167+
label_field='name', geometry_id='road_minor',
168+
color=(0.55, 0.55, 0.55),
169+
densify=False,
170+
merge=True,
171+
mesh_cache=Path(__file__).parent / "capetown_roads_minor_mesh.npz",
172+
)
173+
print(f"Placed {info['geometries']} minor road geometries")
174+
175+
except ImportError as e:
176+
print(f"Skipping roads: {e}")
177+
178+
# --- OpenStreetMap water features ---------------------------------------
179+
try:
180+
water_cache = Path(__file__).parent / "capetown_water.geojson"
181+
water_data = fetch_water(
182+
bounds=BOUNDS,
183+
water_type='all',
184+
cache_path=water_cache,
185+
)
186+
187+
major_features = []
188+
minor_features = []
189+
body_features = []
190+
for f in water_data.get('features', []):
191+
ww = (f.get('properties') or {}).get('waterway', '')
192+
nat = (f.get('properties') or {}).get('natural', '')
193+
if ww in _MAJOR_WATER:
194+
major_features.append(f)
195+
elif ww in _MINOR_WATER:
196+
minor_features.append(f)
197+
elif nat == 'water':
198+
body_features.append(f)
199+
else:
200+
minor_features.append(f)
201+
202+
if major_features:
203+
major_fc = {"type": "FeatureCollection", "features": major_features}
204+
with warnings.catch_warnings():
205+
warnings.filterwarnings("ignore", message="place_geojson called before")
206+
major_info = ds.rtx.place_geojson(
207+
major_fc, z='elevation', height=0,
208+
label_field='name', geometry_id='water_major',
209+
color=(0.40, 0.70, 0.95, 2.25),
210+
densify=False,
211+
merge=True,
212+
mesh_cache=Path(__file__).parent / "capetown_water_major_mesh.npz",
213+
)
214+
print(f"Placed {major_info['geometries']} major water features (rivers, canals)")
215+
216+
if minor_features:
217+
minor_fc = {"type": "FeatureCollection", "features": minor_features}
218+
with warnings.catch_warnings():
219+
warnings.filterwarnings("ignore", message="place_geojson called before")
220+
minor_info = ds.rtx.place_geojson(
221+
minor_fc, z='elevation', height=0,
222+
label_field='name', geometry_id='water_minor',
223+
color=(0.50, 0.75, 0.98, 2.25),
224+
densify=False,
225+
merge=True,
226+
mesh_cache=Path(__file__).parent / "capetown_water_minor_mesh.npz",
227+
)
228+
print(f"Placed {minor_info['geometries']} minor water features (streams, drains)")
229+
230+
if body_features:
231+
body_fc = {"type": "FeatureCollection", "features": body_features}
232+
with warnings.catch_warnings():
233+
warnings.filterwarnings("ignore", message="place_geojson called before")
234+
body_info = ds.rtx.place_geojson(
235+
body_fc, z='elevation', height=0.5,
236+
label_field='name', geometry_id='water_body',
237+
color=(0.35, 0.55, 0.88, 2.25),
238+
extrude=True,
239+
merge=True,
240+
mesh_cache=Path(__file__).parent / "capetown_water_body_mesh.npz",
241+
)
242+
print(f"Placed {body_info['geometries']} water bodies (lakes, ponds)")
243+
244+
except ImportError as e:
245+
print(f"Skipping water features: {e}")
246+
247+
# --- NASA FIRMS fire detections (LANDSAT 30 m, last 7 days) -----------
248+
try:
249+
fire_cache = Path(__file__).parent / "capetown_fires.geojson"
250+
fire_data = fetch_firms(
251+
bounds=BOUNDS,
252+
date_span='7d',
253+
cache_path=fire_cache,
254+
crs='EPSG:32734',
255+
)
256+
if fire_data.get('features'):
257+
elev_scale = 0.025
258+
with warnings.catch_warnings():
259+
warnings.filterwarnings("ignore", message="place_geojson called before")
260+
fire_info = ds.rtx.place_geojson(
261+
fire_data, z='elevation', height=20 * elev_scale,
262+
geometry_id='fire',
263+
color=(1.0, 0.25, 0.0, 3.0),
264+
extrude=True,
265+
merge=True,
266+
)
267+
print(f"Placed {fire_info['geometries']} fire detection footprints")
268+
else:
269+
print("No fire detections in the last 7 days")
270+
except Exception as e:
271+
print(f"Skipping fire layer: {e}")
272+
273+
# --- Wind data --------------------------------------------------------
274+
wind = None
275+
try:
276+
from rtxpy import fetch_wind
277+
wind = fetch_wind(BOUNDS, grid_size=15)
278+
except Exception as e:
279+
print(f"Skipping wind: {e}")
280+
281+
print("\nLaunching explore (press G to cycle layers, Shift+W for wind)...\n")
282+
ds.rtx.explore(
283+
z='elevation',
284+
width=2048,
285+
height=1600,
286+
render_scale=0.5,
287+
color_stretch='cbrt',
288+
subsample=4,
289+
wind_data=wind,
290+
)
291+
292+
print("Done")

0 commit comments

Comments
 (0)