Raster Domain Knowledge Oracle — vibespatial-raster

Dtype Preservation

• Rasters store native dtype: uint8, int16, int32, float32, float64
• Operations should preserve input dtype when possible
• Exceptions: operations that produce a different type (e.g., slope -> float64)
• NEVER upcast uint8 -> float64 unless the computation requires it

2. Nodata Handling

Semantics

• Nodata pixels represent "no information" — they are NOT zero or NaN
• Nodata is a specific sentinel value in the raster's native dtype
• Operations must propagate nodata: if any input pixel in an operation is nodata, the output pixel should be nodata

Nodata in Neighborhood Operations (Focal/Stencil)

• If ANY neighbor is nodata, the center pixel result can be:
  • nodata (strict propagation — default for most operations)
  • computed using only valid neighbors (optional, for statistics)
• Edge pixels (outside raster bounds) are treated as nodata or clamped

Nodata in Algebra Operations

• add(a, b): if either pixel is nodata -> result is nodata
• where(condition, a, b): condition can reference nodata pixels
• classify(raster, breaks): nodata pixels remain nodata in output

Nodata Mask Pattern

# On GPU: nullable mask pointer
d_mask = raster.device_nodata_mask()  # CuPy bool array or None
# In kernel: if (nodata_mask != nullptr && nodata_mask[idx]) { output[idx] = nodata_val; return; }

Common Pitfalls

• Integer nodata (e.g., 255 for uint8) is a valid pixel value in other contexts
• Float nodata should be checked with exact equality, not tolerance
• NaN is a special case: for float rasters, NaN may be used as nodata
• When NaN is nodata, use isnan() in kernels, not == comparison

3. Affine Transforms

Definition

The affine transform maps pixel coordinates to world coordinates:

world_x = a * col + b * row + c
world_y = d * col + e * row + f

Where:

• a = pixel width (x-resolution)
• b = rotation (usually 0)
• c = x-coordinate of upper-left corner
• d = rotation (usually 0)
• e = pixel height (negative for north-up, i.e., y decreases downward)
• f = y-coordinate of upper-left corner

Pixel Centers vs Corners

• The affine maps to the upper-left corner of the pixel
• Pixel center is at (col + 0.5, row + 0.5) in pixel space
• This matters for rasterize (point-in-pixel test) and polygonize (contour placement)

Affine Preservation

• Most operations preserve the input affine (same grid)
• Resampling/reprojection changes the affine
• Subsetting changes c and f (origin shifts)
• Resolution changes affect a and e

4. Focal / Neighborhood Operations

Convolution

• Applies a kernel (weight matrix) to each pixel's neighborhood
• Output[i,j] = sum(kernel * neighborhood around i,j)
• Boundary handling: zero-pad (default), reflect, or clamp

Gaussian Blur

• Convolution with a Gaussian kernel
• Kernel size determined by sigma: size = ceil(6 * sigma) | 1 (ensure odd)
• Separable: apply 1D horizontal then 1D vertical for efficiency

Slope and Aspect

• Computed from elevation rasters (DEM)
• Use Horn's method (3x3 neighborhood):

  dz/dx = ((c + 2f + i) - (a + 2d + g)) / (8 * cell_size_x)
  dz/dy = ((g + 2h + i) - (a + 2b + c)) / (8 * cell_size_y)
  slope = arctan(sqrt(dz_dx^2 + dz_dy^2))
  aspect = arctan2(-dz_dy, dz_dx)
  

• Cell size comes from the affine transform (abs(a) for x, abs(e) for y)

Shared Memory Tiling

• Load tile + halo into shared memory
• Halo width = kernel_radius (e.g., 1 for 3x3, 2 for 5x5)
• Tile dims: (TILE_H + 2*halo, TILE_W + 2*halo)
• __syncthreads() after loading, before computation

5. Connected Component Labeling (CCL)

Algorithm: GPU Union-Find

1. Initialize: each foreground pixel gets its own label (= pixel index)
2. Local merge: check 4-connected neighbors, atomicMin to merge labels
3. Pointer jumping: path compression until convergence
4. Relabel: (optional) compact labels to sequential 1..N

Convergence

• Iterate merge + pointer-jump until no pixel changes label
• Check convergence via device-side changed flag (atomicOr)
• Do NOT check convergence by transferring labels to host each iteration

Connectivity

• 4-connected: up, down, left, right (default)
• 8-connected: adds diagonals (use_8_connectivity parameter)
• Background pixels (value == 0 or == nodata) are excluded

Edge Cases

• All-background raster -> all labels = 0 (no components)
• All-foreground raster -> one component (label = 0)
• Single pixel component -> valid, label = pixel_index
• Very large components can take more iterations to converge

6. Morphology Operations

Binary Erosion

• Output pixel = 1 only if ALL structuring element pixels are foreground
• Shrinks foreground regions, removes small features
• Structuring element typically 3x3 cross or 3x3 square

Binary Dilation

• Output pixel = 1 if ANY structuring element pixel is foreground
• Grows foreground regions, fills small gaps
• Dilation of A by B = complement of erosion of complement(A) by B

Opening (Erosion then Dilation)

• Removes small foreground features while preserving shape
• open(A) = dilate(erode(A))

Closing (Dilation then Erosion)

• Removes small holes while preserving shape
• close(A) = erode(dilate(A))

Sieve (Remove Small Components)

• Run CCL, count component sizes, remove components below threshold
• sieve(raster, min_size=10) removes components with < 10 pixels

GPU Implementation

• Erosion/dilation are stencil operations -> shared memory tiling
• Opening/closing = two sequential stencil launches
• No sync needed between launches on same stream

7. Zonal Statistics

Concept

• Given a raster and a zone raster (integer labels), compute statistics per zone: mean, sum, min, max, count, std
• Zones are non-overlapping labeled regions

GPU Approach

1. Sort pixels by zone ID (CCCL radix_sort)
2. Build segment offsets (unique_by_key or histogram)
3. Segmented reduce per statistic (CCCL segmented_reduce)

Edge Cases

• Zone ID = nodata -> exclude from statistics
• Value = nodata -> exclude from statistics
• Empty zone (no valid pixels) -> statistic = NaN or 0
• Very large zones (>1M pixels) -> still efficient with CCCL

8. Rasterize (Vector to Raster)

Algorithm: Point-in-Polygon per Pixel

• For each pixel center, test if it falls inside any input polygon
• Assign the polygon's value to the pixel (or burn a constant)
• GPU: parallelize over pixels, each thread tests against polygon

Pixel Center Convention

• Test point = pixel center = (col + 0.5, row + 0.5) transformed by affine
• A pixel on the boundary of a polygon should be included (follows rasterio convention)

Edge Cases

• Overlapping polygons: last polygon wins (or configurable merge)
• Polygon with holes: winding number test correctly handles holes
• Multipolygon: treat each part independently
• Empty polygon set: all pixels get nodata

9. Polygonize (Raster to Vector)

Algorithm: Marching Squares

• Classify each 2x2 pixel cell into one of 16 cases based on threshold
• Emit line segments for contour boundaries
• Connect segments into polygon rings
• Assign interior value to each polygon

Cell Classification

• 4 corners of 2x2 cell, each above or below threshold
• 4 bits -> 16 cases (including ambiguous saddle points)
• Saddle resolution: use average of 4 corners to disambiguate

Output

• List of (geometry, value) pairs
• Polygons follow the OGC right-hand rule (exterior CCW, holes CW)
• Coordinates in world space (apply affine transform)

Edge Cases

• Uniform raster (all same value): one polygon covering entire extent
• Nodata pixels: treated as "below threshold" or excluded from polygons
• Sub-pixel resolution: vertices at cell edges (col, row) in pixel space

10. Local Algebra Operations

Element-wise Operations

• add(a, b), subtract(a, b), multiply(a, b), divide(a, b)
• apply(raster, func) — apply a Python function element-wise
• where(condition, true_val, false_val) — conditional selection
• classify(raster, breaks, labels) — reclassify by value ranges

Rules

• Both inputs must have same shape (height, width)
• Both inputs should have same affine and CRS (warn if different)
• Nodata propagation: if either input is nodata -> output is nodata
• Division by zero -> nodata (not inf or NaN)
• Output dtype: result of numpy type promotion rules

GPU Path

• Element-wise ops are CuPy-native (cp.add, cp.where, etc.)
• No custom NVRTC needed
• Zero-copy: operate directly on device arrays