eorst 1.0.1

Earth Observation and Remote Sensing Toolkit - library for raster processing pipelines
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

//! I/O methods for RasterDataset.
//!
//! This module contains methods for reading blocks and writing windows.

use crate::core_types::{RasterData, RasterType};
use crate::types::BlockSize;
use crate::rasterdataset::RasterDataset;
use crate::gdal_utils::{open_for_update, read_raster_band, raster_from_size, write_bands_to_file};

use ndarray::{s, Array3};

impl<R> RasterDataset<R>
where
    R: RasterType,
{
    /// Writes a 3D array to a window in the output file.
    #[doc(hidden)]
    pub fn write_window3<T>(
        &self,
        block_index: usize,
        data: Array3<T>,
        out_fn: &std::path::PathBuf,
        na_value: T,
    ) where
        T: RasterType,
    {
        let trimmed = crate::array_ops::trimm_array3(&data, self.blocks[block_index].overlap_size);

        let block_geotransform = self.blocks[block_index].geo_transform;
        let block_size: BlockSize = BlockSize {
            rows: trimmed.shape()[1],
            cols: trimmed.shape()[2],
        };
        let n_bands = data.shape()[0] as isize;
        let epsg_code = self.metadata.epsg_code;
        let na_value = T::from(na_value).unwrap();
        raster_from_size::<T>(
            out_fn,
            block_geotransform,
            epsg_code,
            block_size,
            n_bands,
            na_value,
        );
        let out_ds = open_for_update(out_fn);
        write_bands_to_file(&out_ds, trimmed, (0, 0), (block_size.cols, block_size.rows));
    }

    /// Reads a block of raster data.
    #[doc(hidden)]
    pub fn read_block<T>(&self, block_id: usize) -> RasterData<T>
    where
        T: RasterType,
    {
        let block = &self.blocks[block_id];
        let read_window = block.read_window;

        let rows = read_window.size.rows as usize;
        let cols = read_window.size.cols as usize;
        let data_shape = (
            self.metadata.shape.times,
            self.metadata.shape.layers,
            rows,
            cols,
        );

        let mut data: RasterData<T> = RasterData::zeros(data_shape);

        for layer in &self.metadata.layers {
            let raster_path = std::path::Path::new(&layer.source);
            let window = (read_window.offset.cols, read_window.offset.rows);
            let window_size = (cols, rows);

            let layer_data = read_raster_band::<T>(raster_path, 1, window, window_size);

            let slice = s![
                layer.time_pos,
                layer.layer_pos,
                ..,
                ..
            ];

            data.slice_mut(slice).assign(&layer_data);
        }
        data
    }

    /// Reads a specific layer from a block.
    #[doc(hidden)]
    pub fn read_block_layer_idx<T>(&self, block_id: usize, layer_idx: usize) -> RasterData<T>
    where
        T: RasterType,
    {
        let block = &self.blocks[block_id];
        let read_window = block.read_window;

        let rows = read_window.size.rows as usize;
        let cols = read_window.size.cols as usize;
        let data_shape = (
            self.metadata.shape.times,
            self.metadata.shape.layers,
            rows,
            cols,
        );

        let mut data: RasterData<T> = RasterData::zeros(data_shape);

        let layer = &self.metadata.layers[layer_idx];
        let raster_path = std::path::Path::new(&layer.source);
        let window = (read_window.offset.cols, read_window.offset.rows);
        let window_size = (cols, rows);

        let layer_data = read_raster_band::<T>(raster_path, 1, window, window_size);

        let slice = s![
            layer.time_pos,
            layer.layer_pos,
            ..,
            ..
        ];

        data.slice_mut(slice).assign(&layer_data);

        data
    }
}