You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by pr...@apache.org on 2026/05/14 21:44:41 UTC
(sedona-spatialbench) 01/01: feat: add COG writer with Perlin noise pixel generation
This is an automated email from the ASF dual-hosted git repository.
prantogg pushed a commit to branch pranav/feature/raster-cog-writer
in repository https://gitbox.apache.org/repos/asf/sedona-spatialbench.git
commit 6c5e9a6a1e84eed18357a99a2e1c6d23fb52752e
Author: Pranav Toggi <pr...@gmail.com>
AuthorDate: Thu May 14 14:44:29 2026 -0700
feat: add COG writer with Perlin noise pixel generation
- Add noise module: 2D Perlin noise generator (deterministic, spatially
correlated) for realistic satellite imagery compression ratios
- Add cog module: GDAL-based COG writer using MEM→COG create_copy
workflow, gated behind 'cog-writer' feature flag
- Single-band UInt8, ZSTD compressed, 256x256 internal tiling
- Per-footprint UTM CRS, geotransform from footprint origin
- Seed from (footprint_id, cog_id) ensures deterministic output
- 24 tests passing (5 noise + 3 COG roundtrip/determinism tests)
---
spatialbench-raster/Cargo.toml | 9 ++
spatialbench-raster/src/cog/mod.rs | 246 +++++++++++++++++++++++++++++++++++
spatialbench-raster/src/lib.rs | 9 +-
spatialbench-raster/src/noise/mod.rs | 200 ++++++++++++++++++++++++++++
4 files changed, 461 insertions(+), 3 deletions(-)
diff --git a/spatialbench-raster/Cargo.toml b/spatialbench-raster/Cargo.toml
index 91d0769..2bebf21 100644
--- a/spatialbench-raster/Cargo.toml
+++ b/spatialbench-raster/Cargo.toml
@@ -28,5 +28,14 @@ description = "Raster benchmark data generation for SpatialBench"
keywords = ["spatial", "raster", "cog", "geotiff", "benchmark"]
categories = ["science::geo", "data-structures"]
+[features]
+default = []
+cog-writer = ["dep:gdal"]
+
[dependencies]
spatialbench = { path = "../spatialbench", version = "0.3.0" }
+gdal = { version = "0.19", optional = true, features = ["bindgen"] }
+
+[dev-dependencies]
+gdal = { version = "0.19", features = ["bindgen"] }
+tempfile = "3"
diff --git a/spatialbench-raster/src/cog/mod.rs b/spatialbench-raster/src/cog/mod.rs
new file mode 100644
index 0000000..9f3f100
--- /dev/null
+++ b/spatialbench-raster/src/cog/mod.rs
@@ -0,0 +1,246 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Cloud-Optimized GeoTIFF (COG) writer using GDAL.
+//!
+//! Generates single-band UInt8 COGs with:
+//! - Per-footprint UTM CRS (EPSG:326xx)
+//! - ZSTD compression
+//! - 256×256 internal tiling
+//! - Perlin noise pixel data for realistic compression ratios
+
+use crate::footprint::Footprint;
+use crate::noise::PerlinNoise;
+
+use gdal::raster::RasterCreationOptions;
+use gdal::spatial_ref::SpatialRef;
+use gdal::{DriverManager, GeoTransform};
+
+use std::io;
+use std::path::Path;
+
+/// Configuration for COG generation.
+#[derive(Debug, Clone, Copy)]
+pub struct CogConfig {
+ /// Raster width in pixels.
+ pub width: u32,
+ /// Raster height in pixels.
+ pub height: u32,
+ /// Internal tile size (pixels per side).
+ pub tile_size: u32,
+ /// Pixel resolution in meters (for geotransform).
+ pub resolution: f64,
+ /// Perlin noise frequency (controls spatial detail per tile).
+ pub noise_frequency: f64,
+}
+
+impl Default for CogConfig {
+ fn default() -> Self {
+ Self {
+ width: 1830,
+ height: 1830,
+ resolution: 60.0,
+ tile_size: 256,
+ noise_frequency: 8.0,
+ }
+ }
+}
+
+/// Write a single COG file for a given footprint and scene ID.
+///
+/// The output is a Cloud-Optimized GeoTIFF with:
+/// - Single band, UInt8, ZSTD compressed
+/// - CRS set to the footprint's UTM zone
+/// - Geotransform derived from the footprint's NW corner origin
+/// - Deterministic Perlin noise pixel data seeded from `(footprint.id, cog_id)`
+///
+/// Internally creates a MEM dataset, populates it, then uses the COG driver's
+/// `create_copy` to produce the final file (COG is a create-copy-only driver).
+///
+/// # Errors
+///
+/// Returns `io::Error` if GDAL fails to create the dataset or write raster data.
+pub fn write_cog(
+ config: &CogConfig,
+ footprint: &Footprint,
+ cog_id: u32,
+ output_path: &Path,
+) -> io::Result<()> {
+ // Create in-memory dataset
+ let mem_driver = DriverManager::get_driver_by_name("MEM").map_err(|e| {
+ io::Error::new(
+ io::ErrorKind::NotFound,
+ format!("GDAL MEM driver not available: {e}"),
+ )
+ })?;
+
+ let mut mem_ds = mem_driver
+ .create_with_band_type::<u8, _>("", config.width as usize, config.height as usize, 1)
+ .map_err(|e| io::Error::other(format!("failed to create MEM dataset: {e}")))?;
+
+ // Set CRS to footprint's UTM zone
+ let srs = SpatialRef::from_epsg(footprint.epsg).map_err(|e| {
+ io::Error::other(format!(
+ "failed to create SRS for EPSG:{}: {e}",
+ footprint.epsg
+ ))
+ })?;
+ mem_ds
+ .set_spatial_ref(&srs)
+ .map_err(|e| io::Error::other(format!("failed to set CRS: {e}")))?;
+
+ // Geotransform: origin is NW corner, positive x-res east, negative y-res south
+ let gt: GeoTransform = [
+ footprint.origin.0, // top-left x (easting)
+ config.resolution, // pixel width (meters)
+ 0.0, // rotation
+ footprint.origin.1, // top-left y (northing)
+ 0.0, // rotation
+ -config.resolution, // pixel height (negative = south)
+ ];
+ mem_ds
+ .set_geo_transform(>)
+ .map_err(|e| io::Error::other(format!("failed to set geotransform: {e}")))?;
+
+ // Generate pixel data with Perlin noise
+ let seed = (footprint.id as u64) << 32 | cog_id as u64;
+ let noise = PerlinNoise::new(seed);
+ let pixels = noise.generate_raster(config.width, config.height, config.noise_frequency);
+
+ // Write raster band
+ let mut band = mem_ds
+ .rasterband(1)
+ .map_err(|e| io::Error::other(format!("failed to get raster band: {e}")))?;
+ let mut buffer =
+ gdal::raster::Buffer::new((config.width as usize, config.height as usize), pixels);
+ band.write::<u8>(
+ (0, 0),
+ (config.width as usize, config.height as usize),
+ &mut buffer,
+ )
+ .map_err(|e| io::Error::other(format!("failed to write pixel data: {e}")))?;
+
+ // Create COG via create_copy from the MEM dataset
+ let cog_driver = DriverManager::get_driver_by_name("COG").map_err(|e| {
+ io::Error::new(
+ io::ErrorKind::NotFound,
+ format!("GDAL COG driver not available: {e}"),
+ )
+ })?;
+
+ let cog_options = RasterCreationOptions::from_iter([
+ "COMPRESS=ZSTD",
+ &format!("BLOCKSIZE={}", config.tile_size),
+ "LEVEL=3",
+ ]);
+
+ mem_ds
+ .create_copy(&cog_driver, output_path, &cog_options)
+ .map_err(|e| io::Error::other(format!("failed to create COG: {e}")))?;
+
+ Ok(())
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+ use crate::footprint::{lon_to_utm_zone, lonlat_to_utm, utm_to_lonlat};
+
+ use gdal::Dataset;
+ use tempfile::tempdir;
+
+ fn test_footprint() -> Footprint {
+ let lon = -100.0_f64;
+ let lat = 35.0_f64;
+ let zone = lon_to_utm_zone(lon);
+ let (e, n) = lonlat_to_utm(lon, lat, zone);
+ let step = 109_800.0;
+ let nw = utm_to_lonlat(e, n, zone, true);
+ let se = utm_to_lonlat(e + step, n - step, zone, true);
+ Footprint {
+ id: 0,
+ epsg: 32600 + zone,
+ origin: (e, n),
+ bbox_4326: [
+ nw.0.min(se.0),
+ nw.1.min(se.1),
+ nw.0.max(se.0),
+ nw.1.max(se.1),
+ ],
+ }
+ }
+
+ #[test]
+ fn cog_roundtrip_default_config() {
+ let dir = tempdir().unwrap();
+ let path = dir.path().join("test.tif");
+ let config = CogConfig::default();
+ let fp = test_footprint();
+
+ write_cog(&config, &fp, 0, &path).unwrap();
+
+ // Verify with GDAL
+ let ds = Dataset::open(&path).unwrap();
+ let (w, h) = ds.raster_size();
+ assert_eq!(w, 1830);
+ assert_eq!(h, 1830);
+ assert_eq!(ds.raster_count(), 1);
+
+ // Verify CRS
+ let srs = ds.spatial_ref().unwrap();
+ assert_eq!(srs.auth_code().unwrap(), fp.epsg as i32);
+
+ // Verify geotransform
+ let gt = ds.geo_transform().unwrap();
+ assert!((gt[0] - fp.origin.0).abs() < 1e-6);
+ assert!((gt[3] - fp.origin.1).abs() < 1e-6);
+ assert!((gt[1] - 60.0).abs() < 1e-6);
+ assert!((gt[5] - (-60.0)).abs() < 1e-6);
+ }
+
+ #[test]
+ fn cog_deterministic() {
+ let dir = tempdir().unwrap();
+ let config = CogConfig::default();
+ let fp = test_footprint();
+
+ let path_a = dir.path().join("a.tif");
+ let path_b = dir.path().join("b.tif");
+ write_cog(&config, &fp, 7, &path_a).unwrap();
+ write_cog(&config, &fp, 7, &path_b).unwrap();
+
+ let a = std::fs::read(&path_a).unwrap();
+ let b = std::fs::read(&path_b).unwrap();
+ assert_eq!(a, b);
+ }
+
+ #[test]
+ fn cog_different_seeds_differ() {
+ let dir = tempdir().unwrap();
+ let config = CogConfig::default();
+ let fp = test_footprint();
+
+ let path_a = dir.path().join("a.tif");
+ let path_b = dir.path().join("b.tif");
+ write_cog(&config, &fp, 0, &path_a).unwrap();
+ write_cog(&config, &fp, 1, &path_b).unwrap();
+
+ let a = std::fs::read(&path_a).unwrap();
+ let b = std::fs::read(&path_b).unwrap();
+ assert_ne!(a, b);
+ }
+}
diff --git a/spatialbench-raster/src/lib.rs b/spatialbench-raster/src/lib.rs
index 85ee831..ab00c6a 100644
--- a/spatialbench-raster/src/lib.rs
+++ b/spatialbench-raster/src/lib.rs
@@ -17,10 +17,13 @@
//! Raster benchmark data generation for SpatialBench.
//!
-//! This crate provides the data model, scaling tables, footprint grid, and
-//! topology factoring logic for the SpatialBench raster benchmark. COG writing
-//! (Phase 1b) and STAC catalog generation are in separate crates.
+//! This crate provides the data model, scaling tables, footprint grid,
+//! topology factoring, and COG writing logic for the SpatialBench raster
+//! benchmark.
+#[cfg(feature = "cog-writer")]
+pub mod cog;
pub mod footprint;
+pub mod noise;
pub mod scaling;
pub mod topology;
diff --git a/spatialbench-raster/src/noise/mod.rs b/spatialbench-raster/src/noise/mod.rs
new file mode 100644
index 0000000..af8d5e1
--- /dev/null
+++ b/spatialbench-raster/src/noise/mod.rs
@@ -0,0 +1,200 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements. See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership. The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License. You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied. See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! 2D Perlin noise generator for deterministic raster pixel generation.
+//!
+//! Produces spatially correlated UInt8 pixel values that mimic the spectral
+//! characteristics of real satellite imagery, yielding realistic compression
+//! ratios (~2-4x with DEFLATE/ZSTD).
+//!
+//! The implementation follows Ken Perlin's improved noise (2002) with a
+//! permutation table seeded deterministically from a footprint/COG ID pair.
+
+/// A seeded 2D Perlin noise generator.
+///
+/// Each instance produces a unique but deterministic noise field based on its
+/// seed. The permutation table is shuffled via a simple LCG seeded from the
+/// input, ensuring reproducibility across runs.
+pub struct PerlinNoise {
+ perm: [u8; 512],
+}
+
+impl PerlinNoise {
+ /// Create a new Perlin noise generator with a deterministic seed.
+ ///
+ /// Typical usage: `seed = (footprint_id as u64) << 32 | cog_id as u64`
+ pub fn new(seed: u64) -> Self {
+ let mut perm = [0u8; 512];
+ // Initialize identity permutation
+ for i in 0..256u16 {
+ perm[i as usize] = i as u8;
+ }
+ // Fisher-Yates shuffle with LCG
+ let mut rng = seed;
+ for i in (1..256usize).rev() {
+ rng = lcg_next(rng);
+ let j = (rng >> 16) as usize % (i + 1);
+ perm.swap(i, j);
+ }
+ // Duplicate for wrapping
+ for i in 0..256 {
+ perm[256 + i] = perm[i];
+ }
+ Self { perm }
+ }
+
+ /// Sample noise at (x, y), returning a value in [-1.0, 1.0].
+ #[inline]
+ pub fn sample(&self, x: f64, y: f64) -> f64 {
+ let xi = x.floor() as i32;
+ let yi = y.floor() as i32;
+ let xf = x - x.floor();
+ let yf = y - y.floor();
+
+ let u = fade(xf);
+ let v = fade(yf);
+
+ // Hash corners
+ let xi = (xi & 255) as usize;
+ let yi = (yi & 255) as usize;
+ let aa = self.perm[self.perm[xi] as usize + yi] as usize;
+ let ab = self.perm[self.perm[xi] as usize + yi + 1] as usize;
+ let ba = self.perm[self.perm[xi + 1] as usize + yi] as usize;
+ let bb = self.perm[self.perm[xi + 1] as usize + yi + 1] as usize;
+
+ // Gradient dot products and bilinear interpolation
+ let x1 = lerp(grad(aa, xf, yf), grad(ba, xf - 1.0, yf), u);
+ let x2 = lerp(grad(ab, xf, yf - 1.0), grad(bb, xf - 1.0, yf - 1.0), u);
+ lerp(x1, x2, v)
+ }
+
+ /// Generate a full raster buffer of UInt8 pixel values.
+ ///
+ /// `width` and `height` are pixel dimensions. `frequency` controls the
+ /// spatial scale of the noise (higher = more detail per tile).
+ /// Returns a row-major buffer of `width * height` bytes.
+ pub fn generate_raster(&self, width: u32, height: u32, frequency: f64) -> Vec<u8> {
+ let len = width as usize * height as usize;
+ let mut buf = Vec::with_capacity(len);
+ for row in 0..height {
+ for col in 0..width {
+ let x = col as f64 * frequency / width as f64;
+ let y = row as f64 * frequency / height as f64;
+ // Sample returns [-1, 1], map to [0, 255]
+ let val = (self.sample(x, y) + 1.0) * 0.5;
+ buf.push((val.clamp(0.0, 1.0) * 255.0) as u8);
+ }
+ }
+ buf
+ }
+}
+
+/// Simple LCG for permutation table shuffling (Numerical Recipes constants).
+#[inline]
+fn lcg_next(state: u64) -> u64 {
+ state
+ .wrapping_mul(6_364_136_223_846_793_005)
+ .wrapping_add(1)
+}
+
+/// Fade curve: 6t^5 - 15t^4 + 10t^3 (Perlin's improved noise, 2002).
+#[inline]
+fn fade(t: f64) -> f64 {
+ t * t * t * (t * (t * 6.0 - 15.0) + 10.0)
+}
+
+/// Linear interpolation.
+#[inline]
+fn lerp(a: f64, b: f64, t: f64) -> f64 {
+ a + t * (b - a)
+}
+
+/// Gradient function using hash to select from 4 directions.
+#[inline]
+fn grad(hash: usize, x: f64, y: f64) -> f64 {
+ match hash & 3 {
+ 0 => x + y,
+ 1 => -x + y,
+ 2 => x - y,
+ _ => -x - y,
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::*;
+
+ #[test]
+ fn deterministic_output() {
+ let n1 = PerlinNoise::new(42);
+ let n2 = PerlinNoise::new(42);
+ for i in 0..100 {
+ let x = i as f64 * 0.1;
+ let y = i as f64 * 0.07;
+ assert_eq!(n1.sample(x, y).to_bits(), n2.sample(x, y).to_bits());
+ }
+ }
+
+ #[test]
+ fn different_seeds_differ() {
+ let n1 = PerlinNoise::new(0);
+ let n2 = PerlinNoise::new(1);
+ let mut differ = false;
+ for i in 0..100 {
+ let x = i as f64 * 0.1;
+ if n1.sample(x, 0.5) != n2.sample(x, 0.5) {
+ differ = true;
+ break;
+ }
+ }
+ assert!(differ, "different seeds should produce different noise");
+ }
+
+ #[test]
+ fn output_range() {
+ let n = PerlinNoise::new(123);
+ for i in 0..1000 {
+ let x = i as f64 * 0.037;
+ let y = i as f64 * 0.053;
+ let val = n.sample(x, y);
+ assert!(
+ val >= -1.0 && val <= 1.0,
+ "sample out of range: {val} at ({x}, {y})"
+ );
+ }
+ }
+
+ #[test]
+ fn generate_raster_correct_size() {
+ let n = PerlinNoise::new(7);
+ let buf = n.generate_raster(64, 64, 4.0);
+ assert_eq!(buf.len(), 64 * 64);
+ }
+
+ #[test]
+ fn generate_raster_has_variation() {
+ let n = PerlinNoise::new(99);
+ let buf = n.generate_raster(128, 128, 8.0);
+ let min = *buf.iter().min().unwrap();
+ let max = *buf.iter().max().unwrap();
+ // With Perlin noise at this frequency, we expect meaningful range
+ assert!(
+ max - min > 50,
+ "expected pixel variation, got range [{min}, {max}]"
+ );
+ }
+}