|
| 1 | +"""Generate a GIF from the playground viewshed/hillshade animation. |
| 2 | +
|
| 3 | +This script creates frames from the Crater Lake hiking animation and |
| 4 | +combines them into a GIF suitable for the README. |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import cupy |
| 9 | +import xarray as xr |
| 10 | +from pathlib import Path |
| 11 | +from PIL import Image |
| 12 | +import io |
| 13 | + |
| 14 | +from rtxpy import RTX, viewshed, hillshade |
| 15 | + |
| 16 | + |
| 17 | +def load_terrain(): |
| 18 | + """Load Crater Lake terrain data.""" |
| 19 | + import rioxarray as rxr |
| 20 | + |
| 21 | + dem_path = Path(__file__).parent / "crater_lake_national_park.tif" |
| 22 | + |
| 23 | + if not dem_path.exists(): |
| 24 | + raise FileNotFoundError(f"DEM file not found at {dem_path}. Run playground.py first to download it.") |
| 25 | + |
| 26 | + print(f"Loading DEM: {dem_path}") |
| 27 | + terrain = rxr.open_rasterio(str(dem_path), masked=True).squeeze() |
| 28 | + |
| 29 | + # Subsample aggressively for smaller GIF file size |
| 30 | + terrain = terrain[::10, ::10] |
| 31 | + |
| 32 | + # Crop edges to remove invalid border values |
| 33 | + crop = 20 |
| 34 | + terrain = terrain[crop:-crop, crop:-crop] |
| 35 | + |
| 36 | + # Scale down elevation for visualization |
| 37 | + terrain.data = terrain.data * 0.2 |
| 38 | + |
| 39 | + # Ensure contiguous array before GPU transfer |
| 40 | + terrain.data = np.ascontiguousarray(terrain.data) |
| 41 | + |
| 42 | + # Convert to cupy for GPU processing |
| 43 | + terrain.data = cupy.asarray(terrain.data) |
| 44 | + |
| 45 | + print(f"Terrain loaded: {terrain.shape}") |
| 46 | + return terrain |
| 47 | + |
| 48 | + |
| 49 | +def generate_hiking_path(x_coords, y_coords, num_points=360): |
| 50 | + """Generate a hiking path around Crater Lake (roughly circular).""" |
| 51 | + cx = (x_coords.min() + x_coords.max()) / 2 |
| 52 | + cy = (y_coords.min() + y_coords.max()) / 2 |
| 53 | + |
| 54 | + rx = (x_coords.max() - x_coords.min()) * 0.25 |
| 55 | + ry = (y_coords.max() - y_coords.min()) * 0.25 |
| 56 | + |
| 57 | + angles = np.linspace(0, 2 * np.pi, num_points) |
| 58 | + wobble = np.sin(angles * 8) * 0.1 |
| 59 | + |
| 60 | + path_x = cx + (rx + rx * wobble) * np.cos(angles) |
| 61 | + path_y = cy + (ry + ry * wobble) * np.sin(angles) |
| 62 | + |
| 63 | + return path_x, path_y |
| 64 | + |
| 65 | + |
| 66 | +def coords_to_pixel(x, y, x_coords, y_coords): |
| 67 | + """Convert data coordinates to pixel coordinates.""" |
| 68 | + px = np.searchsorted(x_coords, x) |
| 69 | + py = np.searchsorted(-y_coords, -y) |
| 70 | + return int(np.clip(px, 0, len(x_coords) - 1)), int(np.clip(py, 0, len(y_coords) - 1)) |
| 71 | + |
| 72 | + |
| 73 | +def draw_legend(colors, x=10, y=10): |
| 74 | + """Draw a legend in the corner of the frame.""" |
| 75 | + from PIL import Image as PILImage, ImageDraw, ImageFont |
| 76 | + |
| 77 | + H, W = colors.shape[:2] |
| 78 | + |
| 79 | + # Create a small PIL image for drawing text |
| 80 | + legend_w, legend_h = 90, 52 |
| 81 | + legend = PILImage.new('RGBA', (legend_w, legend_h), (0, 0, 0, 180)) |
| 82 | + draw = ImageDraw.Draw(legend) |
| 83 | + |
| 84 | + # Use default font |
| 85 | + try: |
| 86 | + font = ImageFont.truetype("/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf", 10) |
| 87 | + except (OSError, IOError): |
| 88 | + font = ImageFont.load_default() |
| 89 | + |
| 90 | + # Legend entries: color swatch + label |
| 91 | + entries = [ |
| 92 | + ((50, 220, 50), "Visible"), |
| 93 | + ((80, 80, 85), "Not Visible"), |
| 94 | + ((0, 255, 255), "Observer"), |
| 95 | + ] |
| 96 | + |
| 97 | + for i, (color, label) in enumerate(entries): |
| 98 | + row_y = 5 + i * 15 |
| 99 | + # Draw color swatch |
| 100 | + draw.rectangle([5, row_y, 15, row_y + 10], fill=color) |
| 101 | + # Draw label |
| 102 | + draw.text((20, row_y - 1), label, fill=(255, 255, 255), font=font) |
| 103 | + |
| 104 | + # Convert legend to numpy and overlay on frame |
| 105 | + legend_arr = np.array(legend) |
| 106 | + |
| 107 | + # Blend legend onto frame |
| 108 | + for ly in range(legend_h): |
| 109 | + for lx in range(legend_w): |
| 110 | + fy, fx = y + ly, x + lx |
| 111 | + if 0 <= fy < H and 0 <= fx < W: |
| 112 | + alpha = legend_arr[ly, lx, 3] / 255.0 |
| 113 | + colors[fy, fx, :3] = ( |
| 114 | + colors[fy, fx, :3] * (1 - alpha) + legend_arr[ly, lx, :3] * alpha |
| 115 | + ).astype(np.uint8) |
| 116 | + |
| 117 | + |
| 118 | +def draw_observer_marker(colors, px, py, radius=6, glow_radius=18): |
| 119 | + """Draw a glowing teal marker with dark outline at the observer's position.""" |
| 120 | + H, W = colors.shape[:2] |
| 121 | + outline_width = 2 |
| 122 | + |
| 123 | + # Draw outer glow, dark outline, then bright center |
| 124 | + for dy in range(-glow_radius, glow_radius + 1): |
| 125 | + for dx in range(-glow_radius, glow_radius + 1): |
| 126 | + dist_sq = dx*dx + dy*dy |
| 127 | + if dist_sq <= glow_radius * glow_radius: |
| 128 | + ny, nx = py + dy, px + dx |
| 129 | + if 0 <= ny < H and 0 <= nx < W: |
| 130 | + dist = np.sqrt(dist_sq) |
| 131 | + if dist <= radius: |
| 132 | + # Bright cyan/teal center |
| 133 | + colors[ny, nx] = [0, 255, 255, 255] |
| 134 | + elif dist <= radius + outline_width: |
| 135 | + # Dark outline for contrast |
| 136 | + colors[ny, nx] = [0, 40, 40, 255] |
| 137 | + elif dist <= glow_radius: |
| 138 | + # Glow falloff - blend cyan with existing color |
| 139 | + t = (dist - radius - outline_width) / (glow_radius - radius - outline_width) |
| 140 | + glow_strength = (1 - t) ** 1.5 # Slightly softer falloff |
| 141 | + existing = colors[ny, nx, :3].astype(np.float32) |
| 142 | + cyan = np.array([0, 255, 255], dtype=np.float32) |
| 143 | + blended = existing + (cyan - existing) * glow_strength * 0.7 |
| 144 | + colors[ny, nx, :3] = np.clip(blended, 0, 255).astype(np.uint8) |
| 145 | + |
| 146 | + |
| 147 | +def generate_frames(terrain, num_frames=72): |
| 148 | + """Generate animation frames. |
| 149 | +
|
| 150 | + Parameters |
| 151 | + ---------- |
| 152 | + terrain : xarray.DataArray |
| 153 | + The terrain data. |
| 154 | + num_frames : int |
| 155 | + Number of frames to generate. Both the hillshade and hiker will |
| 156 | + complete exactly one full 360° loop in this many frames. |
| 157 | + """ |
| 158 | + H, W = terrain.data.shape |
| 159 | + rtx = RTX() |
| 160 | + |
| 161 | + x_coords = terrain.indexes.get('x').values |
| 162 | + y_coords = terrain.indexes.get('y').values |
| 163 | + |
| 164 | + path_x, path_y = generate_hiking_path(x_coords, y_coords, num_points=360) |
| 165 | + |
| 166 | + frames = [] |
| 167 | + azimuth = 225 |
| 168 | + |
| 169 | + print(f"Generating {num_frames} frames...") |
| 170 | + |
| 171 | + # Calculate rotation per frame for full 360° loop |
| 172 | + azimuth_step = 360 / num_frames |
| 173 | + hiker_step = 360 / num_frames |
| 174 | + |
| 175 | + for frame_idx in range(num_frames): |
| 176 | + path_idx = int((frame_idx * hiker_step) % 360) |
| 177 | + |
| 178 | + vsw = path_x[path_idx] |
| 179 | + vsh = path_y[path_idx] |
| 180 | + azimuth = (225 + frame_idx * azimuth_step) % 360 |
| 181 | + |
| 182 | + # Compute hillshade and viewshed |
| 183 | + hs = hillshade(terrain, |
| 184 | + shadows=True, |
| 185 | + azimuth=azimuth, |
| 186 | + angle_altitude=25, |
| 187 | + rtx=rtx) |
| 188 | + vs = viewshed(terrain, |
| 189 | + x=vsw, |
| 190 | + y=vsh, |
| 191 | + observer_elev=100.0, |
| 192 | + rtx=rtx) |
| 193 | + |
| 194 | + # Convert to numpy arrays |
| 195 | + hs_data = hs.data.get() if hasattr(hs.data, 'get') else hs.data |
| 196 | + vs_data = vs.data.get() if hasattr(vs.data, 'get') else vs.data |
| 197 | + |
| 198 | + # Track NaN and zero pixels before converting - these will be transparent |
| 199 | + transparent_mask = np.isnan(hs_data) | np.isnan(vs_data) | (hs_data == 0) |
| 200 | + |
| 201 | + hs_data = np.nan_to_num(hs_data, nan=0.5) |
| 202 | + gray = np.uint8(np.clip(hs_data * 200, 0, 255)) |
| 203 | + |
| 204 | + # Viewshed returns -1 for invisible, 0-180 for visible (angle) |
| 205 | + visible_mask = vs_data >= 0 |
| 206 | + not_visible_mask = (vs_data < 0) & ~transparent_mask |
| 207 | + |
| 208 | + # Compose the final image with alpha channel (RGBA) |
| 209 | + colors = np.zeros((H, W, 4), dtype=np.uint8) |
| 210 | + colors[:, :, 0] = gray |
| 211 | + colors[:, :, 1] = gray |
| 212 | + colors[:, :, 2] = gray |
| 213 | + colors[:, :, 3] = 255 # Fully opaque by default |
| 214 | + |
| 215 | + # Tint visible areas bright lime green - make it really pop! |
| 216 | + colors[visible_mask, 0] = 50 # Low red |
| 217 | + colors[visible_mask, 1] = np.minimum(255, gray[visible_mask].astype(np.int16) + 120).astype(np.uint8) # Bright green |
| 218 | + colors[visible_mask, 2] = 50 # Low blue |
| 219 | + |
| 220 | + # Tint non-visible areas darker gray |
| 221 | + colors[not_visible_mask, 0] = (colors[not_visible_mask, 0] * 0.5).astype(np.uint8) |
| 222 | + colors[not_visible_mask, 1] = (colors[not_visible_mask, 1] * 0.5).astype(np.uint8) |
| 223 | + colors[not_visible_mask, 2] = (colors[not_visible_mask, 2] * 0.55).astype(np.uint8) |
| 224 | + |
| 225 | + # Make NaN and zero pixels transparent |
| 226 | + colors[transparent_mask, 3] = 0 |
| 227 | + |
| 228 | + # Draw observer marker |
| 229 | + px, py = coords_to_pixel(vsw, vsh, x_coords, y_coords) |
| 230 | + draw_observer_marker(colors, px, py, radius=4) |
| 231 | + |
| 232 | + # Draw legend |
| 233 | + draw_legend(colors, x=10, y=10) |
| 234 | + |
| 235 | + frames.append(Image.fromarray(colors, mode='RGBA')) |
| 236 | + |
| 237 | + if (frame_idx + 1) % 10 == 0: |
| 238 | + print(f" Frame {frame_idx + 1}/{num_frames}") |
| 239 | + |
| 240 | + return frames |
| 241 | + |
| 242 | + |
| 243 | +def create_gif(frames, output_path, fps=12, max_colors=64): |
| 244 | + """Create a GIF from frames. |
| 245 | +
|
| 246 | + Parameters |
| 247 | + ---------- |
| 248 | + frames : list of PIL.Image |
| 249 | + The frames to combine. |
| 250 | + output_path : Path |
| 251 | + Output path for the GIF. |
| 252 | + fps : int |
| 253 | + Frames per second. |
| 254 | + max_colors : int |
| 255 | + Maximum colors in palette for smaller file size. |
| 256 | + """ |
| 257 | + duration = int(1000 / fps) # Duration in milliseconds |
| 258 | + |
| 259 | + # Use a magenta color as the transparency key (unlikely to appear in terrain) |
| 260 | + transparent_color = (255, 0, 255) |
| 261 | + |
| 262 | + # Convert RGBA frames to RGB, replacing transparent pixels with the key color |
| 263 | + print("Converting frames for GIF transparency...") |
| 264 | + rgb_frames = [] |
| 265 | + for frame in frames: |
| 266 | + arr = np.array(frame) |
| 267 | + rgb = arr[:, :, :3].copy() |
| 268 | + alpha = arr[:, :, 3] |
| 269 | + # Set transparent pixels to the key color |
| 270 | + rgb[alpha == 0] = transparent_color |
| 271 | + rgb_frames.append(Image.fromarray(rgb, mode='RGB')) |
| 272 | + |
| 273 | + # Create global palette from sampled frames to avoid flickering |
| 274 | + print(f"Building global palette from {len(rgb_frames)} frames...") |
| 275 | + sample_step = max(1, len(rgb_frames) // 10) |
| 276 | + sampled = [np.array(rgb_frames[i]) for i in range(0, len(rgb_frames), sample_step)] |
| 277 | + combined = np.concatenate([p.reshape(-1, 3) for p in sampled], axis=0) |
| 278 | + |
| 279 | + h, w = np.array(rgb_frames[0]).shape[:2] |
| 280 | + sample_h = int(np.ceil(len(combined) / w)) |
| 281 | + padded = np.zeros((sample_h * w, 3), dtype=np.uint8) |
| 282 | + padded[:len(combined)] = combined |
| 283 | + palette_img = Image.fromarray(padded.reshape(sample_h, w, 3), mode='RGB') |
| 284 | + global_palette = palette_img.quantize(colors=max_colors, method=Image.Quantize.MEDIANCUT) |
| 285 | + |
| 286 | + # Modify palette to reserve index 0 for transparency |
| 287 | + palette_data = list(global_palette.getpalette()) |
| 288 | + palette_data[0:3] = transparent_color # Force index 0 to be transparent color |
| 289 | + global_palette.putpalette(palette_data) |
| 290 | + transparency_index = 0 |
| 291 | + |
| 292 | + # Quantize all frames using the global palette |
| 293 | + print(f"Quantizing frames to {max_colors} colors...") |
| 294 | + quantized_frames = [] |
| 295 | + for frame in rgb_frames: |
| 296 | + q_frame = frame.quantize(palette=global_palette, dither=Image.Dither.FLOYDSTEINBERG) |
| 297 | + quantized_frames.append(q_frame) |
| 298 | + |
| 299 | + print(f"Creating GIF at {output_path}...") |
| 300 | + save_kwargs = { |
| 301 | + 'save_all': True, |
| 302 | + 'append_images': quantized_frames[1:], |
| 303 | + 'duration': duration, |
| 304 | + 'loop': 0, # Loop forever |
| 305 | + 'optimize': True |
| 306 | + } |
| 307 | + if transparency_index is not None: |
| 308 | + save_kwargs['transparency'] = transparency_index |
| 309 | + save_kwargs['disposal'] = 2 # Restore to background |
| 310 | + print(f" Using transparency index: {transparency_index}") |
| 311 | + |
| 312 | + quantized_frames[0].save(output_path, **save_kwargs) |
| 313 | + |
| 314 | + file_size = output_path.stat().st_size / (1024 * 1024) |
| 315 | + print(f"GIF created: {output_path} ({file_size:.1f} MB)") |
| 316 | + |
| 317 | + |
| 318 | +def main(): |
| 319 | + output_path = Path(__file__).parent / "images" / "playground_demo.gif" |
| 320 | + output_path.parent.mkdir(exist_ok=True) |
| 321 | + |
| 322 | + terrain = load_terrain() |
| 323 | + |
| 324 | + # Generate 120 frames - both hillshade and hiker complete exactly one 360° loop |
| 325 | + # At 15fps this gives an 8 second loop that repeats seamlessly |
| 326 | + frames = generate_frames(terrain, num_frames=120) |
| 327 | + |
| 328 | + create_gif(frames, output_path, fps=6, max_colors=128) |
| 329 | + |
| 330 | + print(f"\nDone! GIF saved to: {output_path}") |
| 331 | + |
| 332 | + |
| 333 | +if __name__ == "__main__": |
| 334 | + main() |
0 commit comments