oxigdal-vrt 0.1.7

VRT (Virtual Raster) driver for OxiGDAL - Pure Rust GDAL reimplementation
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
//! VRT reader with lazy evaluation

use crate::band::PixelFunction;
use crate::dataset::VrtDataset;
use crate::error::{Result, VrtError};
use crate::mosaic::MosaicCompositor;
use crate::source::{PixelRect, VrtSource};
use crate::xml::VrtXmlParser;
use lru::LruCache;
use oxigdal_core::buffer::RasterBuffer;
use oxigdal_core::io::FileDataSource;
use oxigdal_core::types::{GeoTransform, NoDataValue, RasterDataType, RasterMetadata};
use oxigdal_geotiff::GeoTiffReader;
use std::num::NonZeroUsize;
use std::path::{Path, PathBuf};
use std::sync::{Arc, Mutex};

/// VRT reader with lazy source loading
pub struct VrtReader {
    /// VRT dataset definition
    dataset: VrtDataset,
    /// Cache of opened source datasets
    source_cache: Arc<Mutex<LruCache<PathBuf, Arc<SourceDataset>>>>,
    /// Mosaic compositor
    compositor: MosaicCompositor,
}

impl VrtReader {
    /// Opens a VRT file
    ///
    /// # Errors
    /// Returns an error if the file cannot be opened or parsed
    pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
        let dataset = VrtXmlParser::parse_file(&path)?;
        Self::from_dataset(dataset)
    }

    /// Creates a reader from a VRT dataset
    ///
    /// # Errors
    /// Returns an error if the dataset is invalid
    pub fn from_dataset(dataset: VrtDataset) -> Result<Self> {
        dataset.validate()?;

        // Create source cache (default 32 open files)
        let cache_size =
            NonZeroUsize::new(32).ok_or_else(|| VrtError::cache_error("Invalid cache size"))?;
        let source_cache = Arc::new(Mutex::new(LruCache::new(cache_size)));

        let compositor = MosaicCompositor::new();

        Ok(Self {
            dataset,
            source_cache,
            compositor,
        })
    }

    /// Gets the raster width
    #[must_use]
    pub fn width(&self) -> u64 {
        self.dataset.raster_x_size
    }

    /// Gets the raster height
    #[must_use]
    pub fn height(&self) -> u64 {
        self.dataset.raster_y_size
    }

    /// Gets the number of bands
    #[must_use]
    pub fn band_count(&self) -> usize {
        self.dataset.band_count()
    }

    /// Gets the GeoTransform
    #[must_use]
    pub fn geo_transform(&self) -> Option<&GeoTransform> {
        self.dataset.geo_transform.as_ref()
    }

    /// Gets the spatial reference system
    #[must_use]
    pub fn srs(&self) -> Option<&str> {
        self.dataset.srs.as_deref()
    }

    /// Gets the block size
    #[must_use]
    pub fn block_size(&self) -> (u32, u32) {
        self.dataset.effective_block_size()
    }

    /// Gets the metadata
    #[must_use]
    pub fn metadata(&self) -> RasterMetadata {
        let (tile_width, tile_height) = self.block_size();
        RasterMetadata {
            width: self.dataset.raster_x_size,
            height: self.dataset.raster_y_size,
            band_count: self.dataset.band_count() as u32,
            data_type: self
                .dataset
                .primary_data_type()
                .unwrap_or(RasterDataType::UInt8),
            geo_transform: self.dataset.geo_transform,
            crs_wkt: self.dataset.srs.clone(),
            nodata: NoDataValue::None,
            color_interpretation: Vec::new(),
            layout: oxigdal_core::types::PixelLayout::Tiled {
                tile_width,
                tile_height,
            },
            driver_metadata: Vec::new(),
            statistics: None,
        }
    }

    /// Reads a band's data for a specific window
    ///
    /// # Errors
    /// Returns an error if reading fails
    pub fn read_window(&self, band: usize, window: PixelRect) -> Result<RasterBuffer> {
        let band_idx = band - 1;
        let vrt_band = self
            .dataset
            .get_band(band_idx)
            .ok_or_else(|| VrtError::band_out_of_range(band, self.dataset.band_count()))?;

        // Get sources that intersect with the window
        let contributing_sources: Vec<&VrtSource> = vrt_band
            .sources
            .iter()
            .filter(|s| s.dst_rect().map(|r| r.intersects(&window)).unwrap_or(false))
            .collect();

        if contributing_sources.is_empty() {
            return Err(VrtError::invalid_window(
                "No sources contribute to this window",
            ));
        }

        // Create output buffer
        let data_size = (window.x_size * window.y_size) as usize * vrt_band.data_type.size_bytes();
        let mut data = vec![0u8; data_size];

        // If pixel function is present, read all sources separately and apply function
        if let Some(ref pixel_func) = vrt_band.pixel_function {
            self.apply_pixel_function(
                &contributing_sources,
                &window,
                vrt_band.data_type,
                vrt_band.nodata,
                pixel_func,
                &mut data,
            )?;
        } else {
            // Composite data from all contributing sources (no pixel function)
            for source in &contributing_sources {
                self.read_source_contribution(source, &window, vrt_band.data_type, &mut data)?;
            }
        }

        RasterBuffer::new(
            data,
            window.x_size,
            window.y_size,
            vrt_band.data_type,
            vrt_band.nodata,
        )
        .map_err(|e| e.into())
    }

    /// Reads a full band
    ///
    /// # Errors
    /// Returns an error if reading fails
    pub fn read_band(&self, band: usize) -> Result<RasterBuffer> {
        let window = PixelRect::new(0, 0, self.width(), self.height());
        self.read_window(band, window)
    }

    /// Reads a source's contribution to a window
    fn read_source_contribution(
        &self,
        source: &VrtSource,
        dst_window: &PixelRect,
        data_type: RasterDataType,
        output: &mut [u8],
    ) -> Result<()> {
        let source_dst_rect = source
            .dst_rect()
            .ok_or_else(|| VrtError::invalid_source("Source has no destination rectangle"))?;

        // Calculate intersection between source and requested window
        let intersection = source_dst_rect
            .intersect(dst_window)
            .ok_or_else(|| VrtError::invalid_window("Source does not intersect window"))?;

        // Open source dataset
        let dataset = self.open_source(source)?;

        // Calculate source rectangle
        let src_window = source
            .window
            .as_ref()
            .ok_or_else(|| VrtError::invalid_source("Source has no window configuration"))?;

        // Calculate offset within source
        let src_x_off = src_window.src_rect.x_off + (intersection.x_off - source_dst_rect.x_off);
        let src_y_off = src_window.src_rect.y_off + (intersection.y_off - source_dst_rect.y_off);

        let src_rect = PixelRect::new(
            src_x_off,
            src_y_off,
            intersection.x_size,
            intersection.y_size,
        );

        // Read from source
        let source_data = dataset.read_window(source.source_band, src_rect)?;

        // Copy to output buffer at correct position
        let dst_x_off = intersection.x_off - dst_window.x_off;
        let dst_y_off = intersection.y_off - dst_window.y_off;

        let params = crate::mosaic::CompositeParams::new(
            dst_x_off,
            dst_y_off,
            intersection.x_size,
            intersection.y_size,
            dst_window.x_size,
            data_type,
        );
        self.compositor
            .composite(source_data.as_bytes(), output, &params)?;

        Ok(())
    }

    /// Applies pixel function to source data
    fn apply_pixel_function(
        &self,
        sources: &[&VrtSource],
        window: &PixelRect,
        data_type: RasterDataType,
        nodata: NoDataValue,
        pixel_func: &PixelFunction,
        output: &mut [u8],
    ) -> Result<()> {
        let pixel_count = (window.x_size * window.y_size) as usize;
        let _bytes_per_pixel = data_type.size_bytes();

        // Read all source bands
        let mut source_buffers = Vec::new();
        for source in sources {
            let source_dst_rect = source
                .dst_rect()
                .ok_or_else(|| VrtError::invalid_source("Source has no destination rectangle"))?;

            let intersection = source_dst_rect
                .intersect(window)
                .ok_or_else(|| VrtError::invalid_window("Source does not intersect window"))?;

            let dataset = self.open_source(source)?;

            let src_window = source
                .window
                .as_ref()
                .ok_or_else(|| VrtError::invalid_source("Source has no window configuration"))?;

            let src_x_off =
                src_window.src_rect.x_off + (intersection.x_off - source_dst_rect.x_off);
            let src_y_off =
                src_window.src_rect.y_off + (intersection.y_off - source_dst_rect.y_off);

            let src_rect = PixelRect::new(
                src_x_off,
                src_y_off,
                intersection.x_size,
                intersection.y_size,
            );

            let source_data = dataset.read_window(source.source_band, src_rect)?;
            source_buffers.push((source_data, intersection));
        }

        // Apply pixel function to each pixel
        for pixel_idx in 0..pixel_count {
            let y = pixel_idx as u64 / window.x_size;
            let x = pixel_idx as u64 % window.x_size;
            let global_x = window.x_off + x;
            let global_y = window.y_off + y;

            // Collect values from all sources for this pixel
            let mut values = Vec::new();
            for (source_buffer, intersection) in &source_buffers {
                if global_x >= intersection.x_off
                    && global_x < intersection.x_off + intersection.x_size
                    && global_y >= intersection.y_off
                    && global_y < intersection.y_off + intersection.y_size
                {
                    let local_x = global_x - intersection.x_off;
                    let local_y = global_y - intersection.y_off;
                    let local_idx = (local_y * intersection.x_size + local_x) as usize;

                    // Read value from source buffer
                    let value = self.read_pixel_value(
                        source_buffer.as_bytes(),
                        local_idx,
                        data_type,
                        nodata,
                    )?;
                    values.push(value);
                } else {
                    values.push(None);
                }
            }

            // Apply pixel function
            let result = pixel_func.apply(&values)?;

            // Write result to output
            self.write_pixel_value(output, pixel_idx, result, data_type)?;
        }

        Ok(())
    }

    /// Reads a single pixel value from a buffer
    fn read_pixel_value(
        &self,
        buffer: &[u8],
        pixel_idx: usize,
        data_type: RasterDataType,
        nodata: NoDataValue,
    ) -> Result<Option<f64>> {
        let bytes_per_pixel = data_type.size_bytes();
        let offset = pixel_idx * bytes_per_pixel;

        if offset + bytes_per_pixel > buffer.len() {
            return Ok(None);
        }

        let value = match data_type {
            RasterDataType::UInt8 => buffer[offset] as f64,
            RasterDataType::Int8 => buffer[offset] as i8 as f64,
            RasterDataType::UInt16 => {
                let val = u16::from_le_bytes([buffer[offset], buffer[offset + 1]]);
                val as f64
            }
            RasterDataType::Int16 => {
                let val = i16::from_le_bytes([buffer[offset], buffer[offset + 1]]);
                val as f64
            }
            RasterDataType::UInt32 => {
                let val = u32::from_le_bytes([
                    buffer[offset],
                    buffer[offset + 1],
                    buffer[offset + 2],
                    buffer[offset + 3],
                ]);
                val as f64
            }
            RasterDataType::Int32 => {
                let val = i32::from_le_bytes([
                    buffer[offset],
                    buffer[offset + 1],
                    buffer[offset + 2],
                    buffer[offset + 3],
                ]);
                val as f64
            }
            RasterDataType::Float32 => {
                let val = f32::from_le_bytes([
                    buffer[offset],
                    buffer[offset + 1],
                    buffer[offset + 2],
                    buffer[offset + 3],
                ]);
                val as f64
            }
            RasterDataType::Float64 => f64::from_le_bytes([
                buffer[offset],
                buffer[offset + 1],
                buffer[offset + 2],
                buffer[offset + 3],
                buffer[offset + 4],
                buffer[offset + 5],
                buffer[offset + 6],
                buffer[offset + 7],
            ]),
            _ => return Err(VrtError::invalid_source("Unsupported data type")),
        };

        // Check for NoData
        let is_nodata = match nodata {
            NoDataValue::None => false,
            NoDataValue::Integer(nd) => (value - nd as f64).abs() < f64::EPSILON,
            NoDataValue::Float(nd) => (value - nd).abs() < f64::EPSILON,
        };

        if is_nodata { Ok(None) } else { Ok(Some(value)) }
    }

    /// Writes a single pixel value to a buffer
    fn write_pixel_value(
        &self,
        buffer: &mut [u8],
        pixel_idx: usize,
        value: Option<f64>,
        data_type: RasterDataType,
    ) -> Result<()> {
        let bytes_per_pixel = data_type.size_bytes();
        let offset = pixel_idx * bytes_per_pixel;

        if offset + bytes_per_pixel > buffer.len() {
            return Err(VrtError::invalid_window("Pixel offset out of bounds"));
        }

        let write_val = value.unwrap_or(0.0);

        match data_type {
            RasterDataType::UInt8 => {
                buffer[offset] = write_val.clamp(0.0, 255.0) as u8;
            }
            RasterDataType::Int8 => {
                buffer[offset] = write_val.clamp(-128.0, 127.0) as i8 as u8;
            }
            RasterDataType::UInt16 => {
                let val = write_val.clamp(0.0, 65535.0) as u16;
                buffer[offset..offset + 2].copy_from_slice(&val.to_le_bytes());
            }
            RasterDataType::Int16 => {
                let val = write_val.clamp(-32768.0, 32767.0) as i16;
                buffer[offset..offset + 2].copy_from_slice(&val.to_le_bytes());
            }
            RasterDataType::UInt32 => {
                let val = write_val.clamp(0.0, u32::MAX as f64) as u32;
                buffer[offset..offset + 4].copy_from_slice(&val.to_le_bytes());
            }
            RasterDataType::Int32 => {
                let val = write_val.clamp(i32::MIN as f64, i32::MAX as f64) as i32;
                buffer[offset..offset + 4].copy_from_slice(&val.to_le_bytes());
            }
            RasterDataType::Float32 => {
                let val = write_val as f32;
                buffer[offset..offset + 4].copy_from_slice(&val.to_le_bytes());
            }
            RasterDataType::Float64 => {
                buffer[offset..offset + 8].copy_from_slice(&write_val.to_le_bytes());
            }
            _ => return Err(VrtError::invalid_source("Unsupported data type")),
        }

        Ok(())
    }

    /// Opens a source dataset (with caching)
    fn open_source(&self, source: &VrtSource) -> Result<Arc<SourceDataset>> {
        let path = if let Some(ref vrt_path) = self.dataset.vrt_path {
            source.filename.resolve(vrt_path)?
        } else {
            source.filename.path.clone()
        };

        // Check cache first
        {
            let mut cache = self
                .source_cache
                .lock()
                .map_err(|_| VrtError::cache_error("Failed to lock source cache"))?;

            if let Some(dataset) = cache.get(&path) {
                return Ok(Arc::clone(dataset));
            }
        }

        // Open new dataset
        let dataset = SourceDataset::open(&path)?;
        let arc_dataset = Arc::new(dataset);

        // Add to cache
        {
            let mut cache = self
                .source_cache
                .lock()
                .map_err(|_| VrtError::cache_error("Failed to lock source cache"))?;
            cache.put(path, Arc::clone(&arc_dataset));
        }

        Ok(arc_dataset)
    }

    /// Clears the source cache
    pub fn clear_cache(&mut self) {
        if let Ok(mut cache) = self.source_cache.lock() {
            cache.clear();
        }
    }

    /// Gets the current cache size
    pub fn cache_size(&self) -> usize {
        self.source_cache
            .lock()
            .map(|cache| cache.len())
            .unwrap_or(0)
    }
}

/// Wrapper for source datasets
pub struct SourceDataset {
    /// GeoTIFF reader (for now, only GeoTIFF sources are supported)
    geotiff: Option<GeoTiffReader<FileDataSource>>,
}

impl SourceDataset {
    /// Opens a source dataset
    ///
    /// # Errors
    /// Returns an error if the file cannot be opened
    pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
        // Try to open as GeoTIFF
        match FileDataSource::open(path.as_ref()) {
            Ok(source) => match GeoTiffReader::open(source) {
                Ok(reader) => Ok(Self {
                    geotiff: Some(reader),
                }),
                Err(e) => Err(VrtError::source_error(
                    path.as_ref().display().to_string(),
                    format!("Failed to open as GeoTIFF: {}", e),
                )),
            },
            Err(e) => Err(VrtError::source_error(
                path.as_ref().display().to_string(),
                format!("Failed to open file: {}", e),
            )),
        }
    }

    /// Reads a window from the source dataset
    ///
    /// # Errors
    /// Returns an error if reading fails
    pub fn read_window(&self, band: usize, window: PixelRect) -> Result<RasterBuffer> {
        if let Some(ref geotiff) = self.geotiff {
            // For now, we read the full band and extract the window
            // A more efficient implementation would read only the necessary tiles
            let full_band = geotiff.read_band(0, band - 1).map_err(|e| {
                VrtError::source_error("source", format!("Failed to read band: {}", e))
            })?;

            // Extract window
            let width = geotiff.width() as usize;
            let height = geotiff.height() as usize;
            let data_type = geotiff.data_type().unwrap_or(RasterDataType::UInt8);
            let bytes_per_pixel = data_type.size_bytes();

            let mut window_data = Vec::new();

            for y in 0..window.y_size {
                let src_y = (window.y_off + y) as usize;
                if src_y >= height {
                    break;
                }

                let src_offset = (src_y * width + window.x_off as usize) * bytes_per_pixel;
                let copy_width = window.x_size.min((width as u64) - window.x_off) as usize;
                let copy_bytes = copy_width * bytes_per_pixel;

                if src_offset + copy_bytes <= full_band.len() {
                    window_data.extend_from_slice(&full_band[src_offset..src_offset + copy_bytes]);
                }
            }

            RasterBuffer::new(
                window_data,
                window.x_size,
                window.y_size,
                data_type,
                geotiff.nodata(),
            )
            .map_err(|e| e.into())
        } else {
            Err(VrtError::source_error(
                "unknown",
                "Unsupported source format",
            ))
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::band::VrtBand;
    use crate::source::VrtSource;

    #[test]
    fn test_vrt_reader_creation() {
        let mut dataset = VrtDataset::new(512, 512);
        let source = VrtSource::simple("/test.tif", 1);
        let band = VrtBand::simple(1, RasterDataType::UInt8, source);
        dataset.add_band(band);

        let reader = VrtReader::from_dataset(dataset);
        assert!(reader.is_ok());
        let r = reader.expect("Should create reader");
        assert_eq!(r.width(), 512);
        assert_eq!(r.height(), 512);
        assert_eq!(r.band_count(), 1);
    }

    #[test]
    fn test_cache() {
        let mut dataset = VrtDataset::new(512, 512);
        let source = VrtSource::simple("/test.tif", 1);
        let band = VrtBand::simple(1, RasterDataType::UInt8, source);
        dataset.add_band(band);

        let mut reader = VrtReader::from_dataset(dataset).expect("Should create reader");
        assert_eq!(reader.cache_size(), 0);

        reader.clear_cache();
        assert_eq!(reader.cache_size(), 0);
    }
}