eorst/
gdal_utils.rs

1//! GDAL utility functions for raster processing.
2//!
3//! This module provides helper functions for GDAL command-line operations
4//! including warping, translating, creating VRTs, and mosaic building.
5
6use crate::metadata::Extent;
7use crate::types::{BlockSize, GeoTransform, ImageResolution};
8use crate::core_types::RasterType;
9use anyhow::Result;
10use anyhow::Context;
11use gdal::{Dataset, DatasetOptions, DriverManager, GdalOpenFlags};
12use gdal::{raster::GdalDataType, raster::RasterCreationOptions, spatial_ref::{CoordTransform, SpatialRef}};
13use gdal::vector::{FieldValue, LayerAccess};
14use kdam::par_tqdm;
15use kdam::rayon::prelude::*;
16use std::env;
17use std::collections::BTreeMap;
18use std::path::{Path, PathBuf};
19use std::process::Command;
20use uuid::Uuid;
21
22// ─── Reusable helpers extracted from duplicated patterns ───
23
24/// Creates a rayon ThreadPool with the given number of CPUs.
25///
26/// Sets the `RAYON_NUM_THREADS` environment variable and builds the pool.
27/// Replaces the duplicated pattern found in 11 locations across the codebase.
28pub fn create_rayon_pool(n_cpus: usize) -> rayon::ThreadPool {
29    unsafe { std::env::set_var("RAYON_NUM_THREADS", n_cpus.to_string()) };
30    rayon::ThreadPoolBuilder::new()
31        .num_threads(n_cpus)
32        .build()
33        .expect("Failed to create rayon thread pool")
34}
35
36/// Extracts the file stem as a `&str` from a path.
37///
38/// Replaces the duplicated pattern found in 17 locations across the codebase.
39pub fn file_stem_str(path: &Path) -> &str {
40    path.file_stem()
41        .expect("Path has no file stem")
42        .to_str()
43        .expect("File stem is not valid UTF-8")
44}
45
46/// Opens a GDAL dataset with update access.
47///
48/// Replaces the duplicated `DatasetOptions { open_flags: GDAL_OF_UPDATE, .. }` pattern
49/// found in `blocks.rs::write_samples_feature()`, `blocks.rs::write3()`,
50/// and `io.rs::write_window3()`.
51pub fn open_for_update(path: &Path) -> Dataset {
52    let opts = DatasetOptions {
53        open_flags: GdalOpenFlags::GDAL_OF_UPDATE,
54        ..DatasetOptions::default()
55    };
56    Dataset::open_ex(path, opts).expect("Failed to open dataset for update")
57}
58
59/// Writes each band of a 3D array to a GDAL dataset.
60///
61/// Replaces the duplicated band-write loop found in `blocks.rs::write3()`
62/// and `io.rs::write_window3()`.
63pub fn write_bands_to_file<T: RasterType>(
64    out_ds: &Dataset,
65    data: ndarray::ArrayView3<T>,
66    write_offset: (isize, isize),
67    write_size: (usize, usize),
68) {
69    use gdal::raster::Buffer;
70    use ndarray::s;
71    for band in 0..data.shape()[0] {
72        let b = (band + 1) as isize;
73        let mut out_band = out_ds.rasterband(b as usize).expect("Failed to get raster band");
74        let data_vec: Vec<T> = data.slice(s![band, .., ..]).into_iter().copied().collect();
75        let mut data_buffer = Buffer::new(write_size, data_vec);
76        out_band
77            .write(write_offset, write_size, &mut data_buffer)
78            .expect("Failed to write band");
79    }
80}
81
82/// Runs a GDAL command-line tool with the given arguments.
83///
84/// Replaces the duplicated `.spawn().expect().wait().expect()` pattern found in 3+ locations.
85pub fn run_gdal_command(argv: &[&str]) {
86    Command::new(argv[0])
87        .args(&argv[1..])
88        .spawn()
89        .expect("failed to start gdal command")
90        .wait()
91        .expect("failed to wait for gdal command");
92}
93
94/// Reads a raster band window into a 2D array.
95///
96/// Replaces the duplicated `Dataset::open` → `rasterband` → `read_as` → `to_array` pattern
97/// found in `io.rs::read_block()`, `io.rs::read_block_layer_idx()`, and
98/// `processing.rs::read_block_hybrid()`.
99pub fn read_raster_band<T: RasterType>(
100    raster_path: &Path,
101    band_index: usize,
102    offset: (isize, isize),
103    window_size: (usize, usize),
104) -> ndarray::Array2<T> {
105    let ds = Dataset::open(raster_path).expect(&format!("Unable to open {:?}", raster_path));
106    let raster_band = ds.rasterband(band_index).expect("Failed to get raster band");
107    let array_size = window_size;
108    let e_resample_alg = None;
109    raster_band
110        .read_as::<T>(offset, window_size, array_size, e_resample_alg)
111        .expect("Failed to read raster band")
112        .to_array()
113        .expect("Failed to convert to array")
114}
115
116/// Assembles processed block files into a final output via mosaic + translate + cleanup.
117///
118/// This replaces the duplicated 4-line pattern (mosaic, translate, remove VRT, remove blocks)
119/// found in 9 locations across the codebase.
120pub fn mosaic_translate_cleanup(
121    collected: &[PathBuf],
122    tmp_file: &Path,
123    out_file: &Path,
124    epsg_code: u32,
125) {
126    mosaic(collected, tmp_file, epsg_code, None, None).expect("Could not mosaic to vrt");
127    translate(tmp_file, out_file).expect("Could not translate to geotiff");
128    std::fs::remove_file(tmp_file).expect("Unable to remove the temporary file");
129    collected
130        .iter()
131        .for_each(|f| std::fs::remove_file(f).expect("Unable to remove file"));
132}
133
134/// Assembles time-step block files into a final output via mosaic + translate + cleanup.
135///
136/// Variant for `Vec<Vec<PathBuf>>` where each inner vec corresponds to a time step.
137pub fn mosaic_translate_cleanup_time_steps(
138    collected: &[Vec<PathBuf>],
139    out_file: &Path,
140    epsg_code: u32,
141    n_times: usize,
142) {
143    par_tqdm!((0..n_times).into_par_iter()).for_each(|time_index| {
144        let mut block_fns = Vec::new();
145        for block in collected.iter() {
146            let block_fn = block[time_index].to_owned();
147            block_fns.push(block_fn);
148        }
149        let tmp_vrt = PathBuf::from(create_temp_file("vrt"));
150        let file_stem = file_stem_str(out_file);
151        let time_fn = out_file.with_file_name(format!("{}_{}.tif", file_stem, time_index));
152        mosaic(&block_fns, &tmp_vrt, epsg_code, None, None).expect("Could not mosaic to vrt");
153        translate(&tmp_vrt, &time_fn).expect("Could not translate to geotiff");
154        std::fs::remove_file(tmp_vrt).expect("Unable to remove the temporary file");
155        block_fns
156            .iter()
157            .for_each(|f| std::fs::remove_file(f).expect("Unable to remove file"));
158    });
159}
160
161/// Creates a temporary file with the given extension.
162///
163/// Uses the `TMP_DIR` environment variable if set, otherwise falls back to `/tmp`.
164/// File names are prefixed with `eorst_` for easy identification.
165pub fn create_temp_file(ext: &str) -> String {
166    let dir = env::var("TMP_DIR").unwrap_or("/tmp".to_string());
167    let dir = Path::new(&dir);
168    let file_name = format!("eorst_{}.{}", Uuid::new_v4().simple(), ext);
169    let file_name = dir.join(file_name);
170    file_name.into_os_string().into_string().unwrap()
171}
172
173/// Warps a raster source to a different EPSG projection.
174pub(crate) fn warp(
175    source: PathBuf,
176    target_resolution: Option<ImageResolution>,
177    target_epsg: u32,
178) -> PathBuf {
179    let new_source = create_temp_file("vrt");
180    let mut args: Vec<String> = vec!["gdalwarp".to_string(), "-q".to_string()];
181    if let Some(tr) = target_resolution {
182        args.extend([
183            "-tr".to_string(),
184            format!("{}", tr.x),
185            format!("{}", tr.y),
186            "-tap".to_string(),
187        ]);
188    }
189    args.extend([
190        "-t_srs".to_string(),
191        format!("EPSG:{}", target_epsg),
192        source.to_string_lossy().to_string(),
193        new_source.clone(),
194    ]);
195    let argv: Vec<&str> = args.iter().map(|s| s.as_str()).collect();
196    run_gdal_command(&argv);
197    PathBuf::from(new_source)
198}
199
200/// Warps a raster source to a different extent and resolution.
201pub(crate) fn warp_with_te_tr(
202    source: PathBuf,
203    target_te: &Extent,
204    target_resolution: ImageResolution,
205) -> PathBuf {
206    let new_source = create_temp_file("vrt");
207    let xmin_s = format!("{}", (target_te.xmin * 100.).round() / 100.);
208    let ymin_s = format!("{}", (target_te.ymin * 100.).round() / 100.);
209    let xmax_s = format!("{}", (target_te.xmax * 100.).round() / 100.);
210    let ymax_s = format!("{}", (target_te.ymax * 100.).round() / 100.);
211    let trx_s = format!("{}", (target_resolution.x * 100.).round() / 100.);
212    let try_s = format!("{}", (target_resolution.y * 100.).round() / 100.);
213    let source_s = source.to_string_lossy();
214    let argv = vec![
215        "gdalwarp", "-q", "-te", &xmin_s, &ymin_s, &xmax_s, &ymax_s,
216        "-tr", &trx_s, &try_s, "-r", "nearest",
217        &source_s, &new_source,
218    ];
219    run_gdal_command(&argv);
220    PathBuf::from(new_source)
221}
222
223/// Changes the resolution of a raster source.
224pub(crate) fn change_res(source: PathBuf, target_resolution: ImageResolution) -> PathBuf {
225    let new_source = create_temp_file("vrt");
226    let trx = format!("{}", target_resolution.x);
227    let try_ = format!("{}", target_resolution.y);
228    let source_s = source.to_string_lossy();
229    let argv = vec![
230        "gdalwarp", "-q", "-tr", &trx, &try_, &source_s, &new_source,
231    ];
232    run_gdal_command(&argv);
233    PathBuf::from(new_source)
234}
235
236/// Extracts a single band from a raster source into a VRT.
237///
238/// Uses GDAL's VRT driver in-process instead of spawning a `gdal_translate`
239/// subprocess, which reduces overhead from ~100ms to ~0.1ms per band.
240pub(crate) fn extract_band(source: &Path, band: usize) -> PathBuf {
241    // Open source to extract metadata
242    let src_ds = match Dataset::open(source) {
243        Ok(ds) => ds,
244        Err(_) => {
245            // Fallback to gdal_translate if we can't open in-process
246            let new_source = create_temp_file("vrt");
247            let argv = &[
248                "gdal_translate",
249                "-b",
250                &format!("{}", band),
251                "-q",
252                source.to_str().unwrap(),
253                &new_source,
254            ];
255            run_gdal_command(argv);
256            return PathBuf::from(new_source);
257        }
258    };
259
260    let (xsize, ysize) = src_ds.raster_size();
261    let gt = src_ds.geo_transform().ok();
262    let srs_wkt = src_ds.spatial_ref().ok().map(|s| s.to_wkt().ok()).flatten();
263
264    let band_idx = band;
265    let band_meta = match src_ds.rasterband(band_idx) {
266        Ok(b) => {
267            let dtype = b.band_type();
268            let no_data = b.no_data_value();
269            (Some(dtype), no_data)
270        }
271        Err(_) => (None, None),
272    };
273
274    // Build VRT XML manually — no subprocess needed
275    let wkt_part = match &srs_wkt {
276        Some(wkt) => format!("  <SRS>{}</SRS>\n", wkt),
277        None => String::new(),
278    };
279    let gt_part = match &gt {
280        Some(arr) => format!(
281            "  <GeoTransform> {:.15}, {:.15}, {:.15}, {:.15}, {:.15}, {:.15} </GeoTransform>\n",
282            arr[0], arr[1], arr[2], arr[3], arr[4], arr[5]
283        ),
284        None => String::new(),
285    };
286    let nd_part = match band_meta.1 {
287        Some(nd) => format!("<NoDataValue>{}</NoDataValue>\n", nd),
288        None => String::new(),
289    };
290
291    let source_abs = std::fs::canonicalize(source)
292        .unwrap_or_else(|_| source.to_path_buf())
293        .to_string_lossy()
294        .to_string();
295
296    // Use Display impl for data type string (e.g., "UInt16", "Float32")
297    let dtype_str = match band_meta.0 {
298        Some(dt) => format!("{}", dt),
299        None => "Unknown".to_string(),
300    };
301
302    let vrt_xml = format!(
303        r#"<VRTDataset rasterXSize="{}" rasterYSize="{}">
304{wkt_part}{gt_part}  <VRTRasterBand dataType="{dtype}" band="1">
305    <SimpleSource>
306      <SourceFilename relativeToVRT="0">{source}</SourceFilename>
307      <SourceBand>{band}</SourceBand>
308      {nd_part}    </SimpleSource>
309  </VRTRasterBand>
310</VRTDataset>"#,
311        xsize, ysize,
312        wkt_part = wkt_part,
313        gt_part = gt_part,
314        dtype = dtype_str,
315        source = source_abs,
316        band = band,
317        nd_part = nd_part,
318    );
319
320    let new_source = PathBuf::from(create_temp_file("vrt"));
321    std::fs::write(&new_source, vrt_xml).expect("Failed to write VRT XML");
322    PathBuf::from(new_source)
323}
324
325/// Creates an empty raster file with the given size and properties.
326pub fn raster_from_size<T>(
327    file_name: &Path,
328    geo_transform: GeoTransform,
329    epsg_code: u32,
330    block_size: BlockSize,
331    n_bands: isize,
332    na_value: T,
333) where
334    T: RasterType,
335{
336    let parent_path = file_name.parent().unwrap();
337    std::fs::create_dir(parent_path).unwrap_or(());
338
339    let size_x = block_size.cols;
340    let size_y = block_size.rows;
341    let driver = DriverManager::get_driver_by_name("GTIFF").unwrap();
342    let options =
343        RasterCreationOptions::from_iter(["COMPRESS=LZW", "BLOCKXSIZE=512", "BLOCKYSIZE=512"]);
344
345    let mut dataset = driver
346        .create_with_band_type_with_options::<T, _>(
347            file_name,
348            size_x,
349            size_y,
350            n_bands as usize,
351            &options,
352        )
353        .unwrap();
354    dataset
355        .set_geo_transform(&geo_transform.to_array())
356        .unwrap();
357    let srs = SpatialRef::from_epsg(epsg_code).unwrap();
358    dataset.set_spatial_ref(&srs).unwrap();
359    for band_index in 1..n_bands + 1 {
360        let mut raster_band = dataset.rasterband(band_index as usize).unwrap();
361        let no_data_f64 = na_value
362            .to_f64()
363            .expect("Failed to convert no_data value to f64");
364        raster_band.set_no_data_value(Some(no_data_f64)).unwrap();
365    }
366}
367
368/// Creates a mosaic from multiple raster files.
369///
370/// If `extent` is provided, each input is warped to that extent via `-te`.
371/// If `resolution` is provided, each input is resampled to that resolution via `-tr`.
372pub fn mosaic(
373    collected: &[PathBuf],
374    tmp_file: &Path,
375    epsg_code: u32,
376    extent: Option<Extent>,
377    resolution: Option<f64>,
378) -> Result<()> {
379    let collected_reproj: Vec<PathBuf> = par_tqdm!(collected
380        .par_iter())
381        .map(|image| {
382            let new_source = create_temp_file("vrt");
383            let epsg_s = format!("EPSG:{}", epsg_code);
384            let image_s = image.to_string_lossy().to_string();
385            let mut argv: Vec<String> = vec![
386                "gdalwarp".into(),
387                "-q".into(),
388                "-t_srs".into(),
389                epsg_s,
390            ];
391            // Add optional extent
392            if let Some(ref ext) = extent {
393                argv.push("-te".into());
394                argv.push(ext.xmin.to_string());
395                argv.push(ext.ymin.to_string());
396                argv.push(ext.xmax.to_string());
397                argv.push(ext.ymax.to_string());
398            }
399            // Add optional resolution
400            if let Some(res) = resolution {
401                argv.push("-tr".into());
402                argv.push(res.to_string());
403                argv.push(res.to_string());
404            }
405            argv.push(image_s);
406            argv.push(new_source.clone());
407            let argv_refs: Vec<&str> = argv.iter().map(|s| s.as_str()).collect();
408            run_gdal_command(&argv_refs);
409            PathBuf::from(new_source)
410        })
411        .collect();
412
413    let mut argv: Vec<String> = vec![
414        "gdalbuildvrt".to_string(),
415        "-q".to_string(),
416        tmp_file.to_string_lossy().to_string(),
417    ];
418    argv.extend(collected_reproj.iter().map(|p| p.to_string_lossy().to_string()));
419    let argv_refs: Vec<&str> = argv.iter().map(|s| s.as_str()).collect();
420    run_gdal_command(&argv_refs);
421
422    Ok(())
423}
424
425/// Mosaics multiple raster files into a single output WITHOUT deleting inputs.
426/// Unlike `mosaic_translate_cleanup`, this preserves the input files.
427pub fn mosaic_keep_inputs(
428    collected: &[PathBuf],
429    out_file: &Path,
430    epsg_code: u32,
431    extent: Option<Extent>,
432    resolution: Option<f64>,
433) -> Result<()> {
434    let tmp_file = PathBuf::from(create_temp_file("vrt"));
435    mosaic(collected, &tmp_file, epsg_code, extent, resolution)?;
436    translate(&tmp_file, out_file)?;
437    std::fs::remove_file(&tmp_file).ok();
438    Ok(())
439}
440
441/// Translates a raster file to a different format using the specified driver.
442pub fn translate_with_driver(tmp_fn: &Path, image_fn: &Path, driver_name: &str) -> Result<()> {
443    let argv = vec![
444        "gdal_translate",
445        "-q",
446        "-of",
447        driver_name,
448        "-co",
449        "BIGTIFF=YES",
450        "-co",
451        "COMPRESS=DEFLATE",
452        "-co",
453        "NUM_THREADS=16",
454        tmp_fn.to_str().unwrap(),
455        image_fn.to_str().unwrap(),
456    ];
457    run_gdal_command(&argv);
458    Ok(())
459}
460
461/// Translates a raster file to GeoTIFF format.
462pub fn translate(tmp_fn: &Path, image_fn: &Path) -> Result<()> {
463    translate_with_driver(tmp_fn, image_fn, "GTiff").unwrap();
464    Ok(())
465}
466
467/// Translates an existing GeoTIFF to a Cloud Optimized GeoTIFF.
468///
469/// Uses GDAL's COG driver in-process via `Dataset::create_copy()` — equivalent
470/// to rioxarray's `driver="COG"` approach. The COG driver handles:
471/// - IFD reordering (overviews before full-res data)
472/// - Ghost area for HTTP range requests
473/// - Tile leader/trailer bytes for parallel HTTP access
474/// - Automatic overview generation via `OVERVIEWS=AUTO`
475///
476/// This replaces the previous subprocess-based `gdal_translate -of COG` which
477/// loaded the entire file into memory (OOM risk for large files).
478///
479/// # Arguments
480/// * `src` - Source GeoTIFF path
481/// * `dst` - Destination COG path
482/// * `compression` - Compression algorithm, e.g. "LZW", "DEFLATE", "ZSTD"
483/// * `overview_resampling` - Resampling method for auto-generated overviews
484pub fn translate_to_cog(
485    src: &Path,
486    dst: &Path,
487    compression: &str,
488    overview_resampling: &str,
489) -> Result<()> {
490    let cog_driver = DriverManager::get_driver_by_name("COG")
491        .context("COG driver not available (requires GDAL 3.1+)")?;
492
493    let src_ds = Dataset::open(src)
494        .with_context(|| format!("Failed to open source GeoTIFF {:?}", src))?;
495
496    let options = RasterCreationOptions::from_iter([
497        format!("COMPRESS={}", compression),
498        "BIGTIFF=YES".to_string(),
499        "OVERVIEWS=AUTO".to_string(),
500        format!("RESAMPLING={}", overview_resampling),
501        "NUM_THREADS=ALL_CPUS".to_string(),
502    ]);
503
504    let dst_str = dst.to_str()
505        .context("Destination path is not valid UTF-8")?;
506
507    src_ds.create_copy(&cog_driver, dst_str, &options)
508        .with_context(|| format!("Failed to create COG at {:?}", dst))?;
509
510    Ok(())
511}
512
513/// Returns the widest GDAL data type among all raster bands in the dataset.
514#[allow(dead_code)]
515pub(crate) fn get_widest_type(source: &PathBuf) -> GdalDataType {
516    use log::warn;
517
518    let dataset = Dataset::open(source).unwrap();
519
520    let mut widest: Option<GdalDataType> = None;
521
522    for i in 1..=dataset.raster_count() {
523        let band = dataset.rasterband(i).expect("Failed to read band");
524        let dtype = band.band_type();
525
526        if let Some(existing) = widest {
527            if dtype != existing {
528                warn!(
529                    "Band {} has different type ({:?}) than first band ({:?})",
530                    i, dtype, existing
531                );
532            }
533            widest = Some(existing.union(dtype));
534        } else {
535            widest = Some(dtype);
536        }
537    }
538
539    widest.expect("Dataset has no bands")
540}
541
542/// Gets class IDs from a vector dataset.
543#[allow(dead_code)]
544pub fn get_class(
545    dataset_path: &PathBuf,
546    id_column: &str,
547    class_column: &str,
548) -> BTreeMap<i16, i16> {
549    let dataset = Dataset::open(dataset_path).unwrap();
550    let mut layer = dataset.layer(0).unwrap();
551
552    let _fields_defn = layer
553        .defn()
554        .fields()
555        .map(|field| (field.name(), field.field_type(), field.width()))
556        .collect::<Vec<_>>();
557    let mut class: BTreeMap<i16, i16> = BTreeMap::new();
558
559    for feature in layer.features() {
560        let id_column_idx = feature.field_index(id_column).expect("Bad column name");
561        let class_column_idx = feature.field_index(class_column).expect("Bad column name");
562        let id = feature
563            .field(id_column_idx)
564            .unwrap()
565            .unwrap()
566            .into_int()
567            .unwrap();
568        let condition = feature
569            .field(class_column_idx)
570            .unwrap()
571            .unwrap_or(FieldValue::IntegerValue(-1))
572            .into_int()
573            .unwrap();
574        class.insert(id as i16, condition as i16);
575    }
576    class
577}
578
579/// Basic metadata extracted from a single GDAL dataset open.
580///
581/// Consolidates the 6 separate `Dataset::open` calls in `builder.rs`
582/// (`get_resolution`, `get_geotransform`, `get_extent`, `get_epsg_code`,
583/// `get_na_value`, `get_image_size`) into a single dataset open.
584#[derive(Debug, Clone)]
585pub struct BasicRasterInfo {
586    /// GeoTransform (x_ul, x_res, x_rot, y_ul, y_rot, y_res)
587    pub geo_transform: [f64; 6],
588    /// Image size in pixels (cols, rows)
589    pub size: (usize, usize),
590    /// EPSG code (0 if unknown)
591    pub epsg_code: u32,
592    /// No-data value (None if not set)
593    pub no_data: Option<f64>,
594    /// Number of raster bands
595    pub n_bands: usize,
596}
597
598impl BasicRasterInfo {
599    /// Resolution (x, y) derived from geo_transform.
600    pub fn resolution(&self) -> crate::types::ImageResolution {
601        crate::types::ImageResolution {
602            x: self.geo_transform[1],
603            y: self.geo_transform[5],
604        }
605    }
606
607    /// GeoTransform struct derived from geo_transform array.
608    pub fn geo_transform_struct(&self) -> crate::types::GeoTransform {
609        crate::types::GeoTransform {
610            x_ul: self.geo_transform[0],
611            x_res: self.geo_transform[1],
612            x_rot: self.geo_transform[2],
613            y_ul: self.geo_transform[3],
614            y_rot: self.geo_transform[4],
615            y_res: self.geo_transform[5],
616        }
617    }
618
619    /// ImageSize struct derived from size tuple.
620    pub fn image_size(&self) -> crate::types::ImageSize {
621        crate::types::ImageSize {
622            rows: self.size.1,
623            cols: self.size.0,
624        }
625    }
626
627    /// No-data value cast to type T, or T::zero() if not set.
628    pub fn na_value<T: crate::core_types::RasterType>(&self) -> T {
629        self.no_data
630            .and_then(|v| num_traits::NumCast::from(v))
631            .unwrap_or(T::zero())
632    }
633}
634
635/// Opens a raster dataset once and extracts all basic metadata fields.
636///
637/// This replaces the pattern of 6 separate `Dataset::open` calls in builder.rs
638/// with a single open, reducing I/O overhead.
639pub fn read_basic_raster_info(source: &Path) -> BasicRasterInfo {
640    let ds = Dataset::open(source).expect(&format!("Unable to open {:?}", source));
641    let geo_transform = ds.geo_transform().expect("Failed to get geo transform");
642    let size = ds.raster_size();
643    let spatial_ref = ds.spatial_ref().expect("Failed to get spatial ref");
644    let epsg_code = spatial_ref.auth_code().unwrap_or(0);
645
646    // Extract no-data value from first band
647    let mut no_data = None;
648    if size.0 > 0 && size.1 > 0 {
649        if let Ok(band) = ds.rasterband(1) {
650            no_data = band.no_data_value();
651        }
652    }
653
654    BasicRasterInfo {
655        geo_transform,
656        size,
657        epsg_code: epsg_code as u32,
658        no_data,
659        n_bands: ds.raster_count(),
660    }
661}
662
663/// Compute the extent of a single raster file in the target CRS.
664///
665/// Reads the GeoTransform and size, computes the four corners, and
666/// reprojects them to the target EPSG using GDAL's CoordTransform API.
667fn compute_single_raster_extent(source: &Path, target_epsg: u32) -> Result<Extent> {
668    let ds = Dataset::open(source)?;
669    let gt = ds.geo_transform()?;
670    let (cols, rows) = ds.raster_size();
671    let src_srs = ds.spatial_ref()?;
672    let src_epsg = src_srs.auth_code().unwrap_or(0) as u32;
673
674    // Compute four corners from GeoTransform: (x_ul, x_res, _, y_ul, _, y_res)
675    let corners = [
676        (gt[0], gt[3]),                                          // top-left
677        (gt[0] + gt[1] * cols as f64, gt[3]),                   // top-right
678        (gt[0] + gt[1] * cols as f64, gt[3] + gt[5] * rows as f64), // bottom-right
679        (gt[0], gt[3] + gt[5] * rows as f64),                   // bottom-left
680    ];
681
682    if src_epsg == target_epsg {
683        // Same CRS — no transform needed
684        let xs: Vec<f64> = corners.iter().map(|(x, _)| *x).collect();
685        let ys: Vec<f64> = corners.iter().map(|(_, y)| *y).collect();
686        Ok(Extent {
687            xmin: xs.iter().cloned().fold(f64::INFINITY, f64::min),
688            ymin: ys.iter().cloned().fold(f64::INFINITY, f64::min),
689            xmax: xs.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
690            ymax: ys.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
691        })
692    } else {
693        // Transform corners to target CRS
694        let target_srs = SpatialRef::from_epsg(target_epsg)?;
695        let ct = CoordTransform::new(&src_srs, &target_srs)?;
696        let mut xs = corners.iter().map(|(x, _)| *x).collect::<Vec<_>>();
697        let mut ys = corners.iter().map(|(_, y)| *y).collect::<Vec<_>>();
698        let mut zs = vec![0.0; xs.len()];
699        ct.transform_coords(&mut xs, &mut ys, &mut zs)?;
700        Ok(Extent {
701            xmin: xs.iter().cloned().fold(f64::INFINITY, f64::min),
702            ymin: ys.iter().cloned().fold(f64::INFINITY, f64::min),
703            xmax: xs.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
704            ymax: ys.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
705        })
706    }
707}
708
709/// Compute the union extent of multiple raster files in the target CRS.
710///
711/// Each input's four corners are reprojected to the target EPSG using
712/// GDAL's CoordTransform API, then the bounding box of all transformed
713/// corners is returned.
714pub fn compute_raster_union_extent(files: &[PathBuf], target_epsg: u32) -> Result<Extent> {
715    let mut union_extent = compute_single_raster_extent(&files[0], target_epsg)?;
716    for file in &files[1..] {
717        let extent = compute_single_raster_extent(file, target_epsg)?;
718        union_extent = union_extent.union(&extent);
719    }
720    Ok(union_extent)
721}
722
723/// Compute the extent of a vector layer in the target CRS.
724///
725/// Reads the layer extent and reprojects the four corners to the target
726/// EPSG using GDAL's CoordTransform API.
727pub fn compute_vector_extent(vector_path: &Path, target_epsg: u32) -> Result<Extent> {
728    let ds = Dataset::open(vector_path)?;
729    let layer = ds.layer(0)?;
730    let ext = layer.get_extent()?;
731
732    let src_srs = layer.spatial_ref().ok_or_else(|| anyhow::anyhow!("Vector layer has no spatial reference"))?;
733    let src_epsg = src_srs.auth_code().unwrap_or(0) as u32;
734
735    let corners = [
736        (ext.MinX, ext.MinY),
737        (ext.MaxX, ext.MinY),
738        (ext.MaxX, ext.MaxY),
739        (ext.MinX, ext.MaxY),
740    ];
741
742    if src_epsg == target_epsg {
743        Ok(Extent {
744            xmin: ext.MinX,
745            ymin: ext.MinY,
746            xmax: ext.MaxX,
747            ymax: ext.MaxY,
748        })
749    } else {
750        let target_srs = SpatialRef::from_epsg(target_epsg)?;
751        let ct = CoordTransform::new(&src_srs, &target_srs)?;
752        let mut xs = corners.iter().map(|(x, _)| *x).collect::<Vec<_>>();
753        let mut ys = corners.iter().map(|(_, y)| *y).collect::<Vec<_>>();
754        let mut zs = vec![0.0; xs.len()];
755        ct.transform_coords(&mut xs, &mut ys, &mut zs)?;
756        Ok(Extent {
757            xmin: xs.iter().cloned().fold(f64::INFINITY, f64::min),
758            ymin: ys.iter().cloned().fold(f64::INFINITY, f64::min),
759            xmax: xs.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
760            ymax: ys.iter().cloned().fold(f64::NEG_INFINITY, f64::max),
761        })
762    }
763}
eorst/gdal_utils.rs

eorst/
gdal_utils.rs
