geotiff_reader/
lib.rs

1//! Pure-Rust GeoTIFF reader with optional HTTP range-backed remote access.
2//!
3//! Supports:
4//! - **GeoTIFF**: TIFF files with GeoKey metadata (EPSG codes, CRS, tiepoints, pixel scale)
5//! - **COG**: overview discovery plus optional remote open via HTTP range requests
6//!
7//! # Example
8//!
9//! ```no_run
10//! # #[cfg(feature = "local")]
11//! # fn main() -> Result<(), geotiff_reader::Error> {
12//! use geotiff_reader::GeoTiffFile;
13//!
14//! let file = GeoTiffFile::open("dem.tif")?;
15//! println!("EPSG: {:?}", file.epsg());
16//! println!("bounds: {:?}", file.geo_bounds());
17//! println!("size: {}x{}", file.width(), file.height());
18//! # Ok(())
19//! # }
20//! # #[cfg(not(feature = "local"))]
21//! # fn main() {}
22//! ```
23
24pub mod crs;
25pub mod error;
26pub mod geokeys;
27pub mod transform;
28
29#[cfg(feature = "cog")]
30pub mod cog;
31
32pub use error::{Error, Result};
33
34#[cfg(feature = "local")]
35use crs::CrsInfo;
36#[cfg(feature = "local")]
37use geokeys::GeoKeyDirectory;
38#[cfg(feature = "local")]
39use ndarray::ArrayD;
40#[cfg(feature = "local")]
41use std::path::Path;
42#[cfg(feature = "local")]
43use tiff_reader::{OpenOptions as TiffOpenOptions, TagValue, TiffFile, TiffSample};
44#[cfg(feature = "local")]
45use transform::GeoTransform;
46
47#[cfg(feature = "local")]
48use geotiff_core::tags::{
49    TAG_GDAL_NODATA, TAG_GEO_ASCII_PARAMS, TAG_GEO_DOUBLE_PARAMS, TAG_GEO_KEY_DIRECTORY,
50    TAG_MODEL_PIXEL_SCALE, TAG_MODEL_TIEPOINT, TAG_MODEL_TRANSFORMATION, TAG_NEW_SUBFILE_TYPE,
51    TAG_SUBFILE_TYPE,
52};
53
54/// A GeoTIFF file handle with geospatial metadata.
55#[cfg(feature = "local")]
56pub struct GeoTiffFile {
57    tiff: TiffFile,
58    geo_metadata: GeoMetadata,
59    crs: CrsInfo,
60    geokeys: GeoKeyDirectory,
61    transform: Option<GeoTransform>,
62    overview_ifds: Vec<usize>,
63}
64
65#[cfg(feature = "local")]
66pub use tiff_reader::OpenOptions as GeoTiffOpenOptions;
67
68pub use geotiff_core::GeoMetadata;
69
70#[cfg(feature = "local")]
71impl GeoTiffFile {
72    /// Open a GeoTIFF file from disk.
73    pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
74        Self::open_with_options(path, TiffOpenOptions::default())
75    }
76
77    /// Open a GeoTIFF file from disk with explicit TIFF decoder options.
78    pub fn open_with_options<P: AsRef<Path>>(path: P, options: GeoTiffOpenOptions) -> Result<Self> {
79        let tiff = TiffFile::open_with_options(path, options)?;
80        Self::from_tiff(tiff)
81    }
82
83    /// Open a GeoTIFF from an owned byte buffer.
84    pub fn from_bytes(data: Vec<u8>) -> Result<Self> {
85        Self::from_bytes_with_options(data, TiffOpenOptions::default())
86    }
87
88    /// Open a GeoTIFF from bytes with explicit TIFF decoder options.
89    pub fn from_bytes_with_options(data: Vec<u8>, options: GeoTiffOpenOptions) -> Result<Self> {
90        let tiff = TiffFile::from_bytes_with_options(data, options)?;
91        Self::from_tiff(tiff)
92    }
93
94    pub(crate) fn from_tiff(tiff: TiffFile) -> Result<Self> {
95        let ifd = tiff.ifd(0)?;
96        let geokeys = parse_geokey_directory(ifd)?;
97        let crs = CrsInfo::from_geokeys(&geokeys);
98        let epsg = crs.epsg();
99        let tiepoints = parse_tiepoints(ifd);
100        let pixel_scale =
101            parse_fixed_len_double_tag::<3>(ifd.tag(TAG_MODEL_PIXEL_SCALE).map(|tag| &tag.value));
102        let transformation = parse_fixed_len_double_tag::<16>(
103            ifd.tag(TAG_MODEL_TRANSFORMATION).map(|tag| &tag.value),
104        );
105        let transform = transformation
106            .as_ref()
107            .map(GeoTransform::from_transformation_matrix)
108            .or_else(|| {
109                let tiepoint = tiepoints.first()?;
110                let scale = pixel_scale.as_ref()?;
111                Some(GeoTransform::from_tiepoint_and_scale_with_raster_type(
112                    tiepoint,
113                    scale,
114                    crs.raster_type_enum(),
115                ))
116            });
117        let geo_bounds = transform
118            .as_ref()
119            .map(|gt| gt.bounds(ifd.width(), ifd.height()));
120        let overview_ifds = tiff
121            .ifds()
122            .iter()
123            .enumerate()
124            .skip(1)
125            .filter_map(|(index, candidate)| is_overview_ifd(ifd, candidate).then_some(index))
126            .collect();
127
128        let geo_metadata = GeoMetadata {
129            epsg,
130            tiepoints,
131            pixel_scale,
132            transformation,
133            nodata: parse_nodata(ifd),
134            band_count: ifd.samples_per_pixel() as u32,
135            width: ifd.width(),
136            height: ifd.height(),
137            geo_bounds,
138        };
139
140        Ok(Self {
141            tiff,
142            geo_metadata,
143            crs,
144            geokeys,
145            transform,
146            overview_ifds,
147        })
148    }
149
150    /// Returns the underlying TIFF file.
151    pub fn tiff(&self) -> &TiffFile {
152        &self.tiff
153    }
154
155    /// Returns the parsed GeoTIFF metadata.
156    pub fn metadata(&self) -> &GeoMetadata {
157        &self.geo_metadata
158    }
159
160    /// Returns the EPSG code of the coordinate reference system, if present.
161    pub fn epsg(&self) -> Option<u32> {
162        self.geo_metadata.epsg
163    }
164
165    /// Returns the extracted CRS information.
166    pub fn crs(&self) -> &CrsInfo {
167        &self.crs
168    }
169
170    /// Returns the parsed GeoKey directory.
171    pub fn geokeys(&self) -> &GeoKeyDirectory {
172        &self.geokeys
173    }
174
175    /// Returns the affine transform, if present.
176    pub fn transform(&self) -> Option<&GeoTransform> {
177        self.transform.as_ref()
178    }
179
180    /// Returns the geographic bounds as `(min_x, min_y, max_x, max_y)`.
181    pub fn geo_bounds(&self) -> Option<[f64; 4]> {
182        self.geo_metadata.geo_bounds
183    }
184
185    /// Convert a pixel coordinate to map coordinates.
186    pub fn pixel_to_geo(&self, col: f64, row: f64) -> Option<(f64, f64)> {
187        self.transform
188            .map(|transform| transform.pixel_to_geo(col, row))
189    }
190
191    /// Convert map coordinates to pixel coordinates.
192    pub fn geo_to_pixel(&self, x: f64, y: f64) -> Option<(f64, f64)> {
193        self.transform
194            .and_then(|transform| transform.geo_to_pixel(x, y))
195    }
196
197    /// Returns the image width in pixels.
198    pub fn width(&self) -> u32 {
199        self.geo_metadata.width
200    }
201
202    /// Returns the image height in pixels.
203    pub fn height(&self) -> u32 {
204        self.geo_metadata.height
205    }
206
207    /// Returns the number of bands.
208    pub fn band_count(&self) -> u32 {
209        self.geo_metadata.band_count
210    }
211
212    /// Returns the nodata value, if set.
213    pub fn nodata(&self) -> Option<&str> {
214        self.geo_metadata.nodata.as_deref()
215    }
216
217    /// Returns the number of internal overview IFDs.
218    pub fn overview_count(&self) -> usize {
219        self.overview_ifds.len()
220    }
221
222    /// Returns the TIFF IFD index of the requested overview.
223    pub fn overview_ifd_index(&self, overview_index: usize) -> Result<usize> {
224        self.overview_ifds
225            .get(overview_index)
226            .copied()
227            .ok_or(Error::OverviewNotFound(overview_index))
228    }
229
230    /// Decode the base-resolution raster into a typed ndarray.
231    pub fn read_raster<T: TiffSample>(&self) -> Result<ArrayD<T>> {
232        self.tiff.read_image::<T>(0).map_err(Into::into)
233    }
234
235    /// Decode a base-resolution pixel window into a typed ndarray.
236    pub fn read_window<T: TiffSample>(
237        &self,
238        row_off: usize,
239        col_off: usize,
240        rows: usize,
241        cols: usize,
242    ) -> Result<ArrayD<T>> {
243        self.tiff
244            .read_window::<T>(0, row_off, col_off, rows, cols)
245            .map_err(Into::into)
246    }
247
248    /// Decode an overview raster into a typed ndarray.
249    pub fn read_overview<T: TiffSample>(&self, overview_index: usize) -> Result<ArrayD<T>> {
250        let ifd_index = self.overview_ifd_index(overview_index)?;
251        self.tiff.read_image::<T>(ifd_index).map_err(Into::into)
252    }
253
254    /// Decode an overview pixel window into a typed ndarray.
255    pub fn read_overview_window<T: TiffSample>(
256        &self,
257        overview_index: usize,
258        row_off: usize,
259        col_off: usize,
260        rows: usize,
261        cols: usize,
262    ) -> Result<ArrayD<T>> {
263        let ifd_index = self.overview_ifd_index(overview_index)?;
264        self.tiff
265            .read_window::<T>(ifd_index, row_off, col_off, rows, cols)
266            .map_err(Into::into)
267    }
268}
269
270#[cfg(feature = "local")]
271fn is_overview_ifd(base: &tiff_reader::Ifd, candidate: &tiff_reader::Ifd) -> bool {
272    let smaller = candidate.width() < base.width() || candidate.height() < base.height();
273    if !smaller {
274        return false;
275    }
276
277    let same_layout = candidate.samples_per_pixel() == base.samples_per_pixel()
278        && candidate.bits_per_sample() == base.bits_per_sample()
279        && candidate.sample_format() == base.sample_format()
280        && candidate.photometric_interpretation() == base.photometric_interpretation();
281    if !same_layout {
282        return false;
283    }
284
285    candidate
286        .tag(TAG_NEW_SUBFILE_TYPE)
287        .and_then(|tag| tag.value.as_u64())
288        .map(|flags| flags & 0x1 != 0)
289        .or_else(|| {
290            candidate
291                .tag(TAG_SUBFILE_TYPE)
292                .and_then(|tag| tag.value.as_u16())
293                .map(|value| value == 2)
294        })
295        .unwrap_or(true)
296}
297
298#[cfg(feature = "local")]
299fn parse_geokey_directory(ifd: &tiff_reader::Ifd) -> Result<GeoKeyDirectory> {
300    let directory = ifd
301        .tag(TAG_GEO_KEY_DIRECTORY)
302        .and_then(|tag| match &tag.value {
303            TagValue::Short(values) => Some(values.as_slice()),
304            _ => None,
305        })
306        .ok_or(Error::NotGeoTiff)?;
307    let double_params = ifd
308        .tag(TAG_GEO_DOUBLE_PARAMS)
309        .and_then(|tag| tag.value.as_f64_vec())
310        .unwrap_or_default();
311    let ascii_params = ifd
312        .tag(TAG_GEO_ASCII_PARAMS)
313        .and_then(|tag| tag.value.as_str())
314        .unwrap_or("");
315    GeoKeyDirectory::parse(directory, &double_params, ascii_params)
316        .ok_or(Error::InvalidGeoKeyDirectory)
317}
318
319#[cfg(feature = "local")]
320fn parse_fixed_len_double_tag<const N: usize>(value: Option<&TagValue>) -> Option<[f64; N]> {
321    let values = value.and_then(TagValue::as_f64_vec)?;
322    if values.len() < N {
323        return None;
324    }
325    let mut out = [0.0; N];
326    out.copy_from_slice(&values[..N]);
327    Some(out)
328}
329
330#[cfg(feature = "local")]
331fn parse_tiepoints(ifd: &tiff_reader::Ifd) -> Vec<[f64; 6]> {
332    let values = ifd
333        .tag(TAG_MODEL_TIEPOINT)
334        .and_then(|tag| tag.value.as_f64_vec())
335        .unwrap_or_default();
336    values
337        .chunks_exact(6)
338        .map(|chunk| [chunk[0], chunk[1], chunk[2], chunk[3], chunk[4], chunk[5]])
339        .collect()
340}
341
342#[cfg(feature = "local")]
343fn parse_nodata(ifd: &tiff_reader::Ifd) -> Option<String> {
344    ifd.tag(TAG_GDAL_NODATA)
345        .and_then(|tag| tag.value.as_str())
346        .map(ToOwned::to_owned)
347}
348
349#[cfg(test)]
350#[cfg(feature = "local")]
351mod tests {
352    use super::GeoTiffFile;
353
354    #[derive(Clone)]
355    struct TestIfdSpec {
356        entries: Vec<(u16, u16, u32, Vec<u8>)>,
357        image_data: Vec<u8>,
358    }
359
360    fn le_u16(value: u16) -> [u8; 2] {
361        value.to_le_bytes()
362    }
363
364    fn le_u32(value: u32) -> [u8; 4] {
365        value.to_le_bytes()
366    }
367
368    fn le_f64(value: f64) -> [u8; 8] {
369        value.to_le_bytes()
370    }
371
372    fn build_classic_tiff(ifds: &[TestIfdSpec]) -> Vec<u8> {
373        let mut ifd_offsets = Vec::with_capacity(ifds.len());
374        let mut cursor = 8usize;
375        for ifd in ifds {
376            ifd_offsets.push(cursor as u32);
377            let deferred_len: usize = ifd
378                .entries
379                .iter()
380                .filter(|(tag, _, _, value)| *tag != 273 && value.len() > 4)
381                .map(|(_, _, _, value)| value.len())
382                .sum();
383            cursor += 2 + ifd.entries.len() * 12 + 4 + ifd.image_data.len() + deferred_len;
384        }
385
386        let mut bytes = Vec::with_capacity(cursor);
387        bytes.extend_from_slice(b"II");
388        bytes.extend_from_slice(&le_u16(42));
389        bytes.extend_from_slice(&le_u32(ifd_offsets.first().copied().unwrap_or(0)));
390
391        for (ifd_index, ifd) in ifds.iter().enumerate() {
392            let ifd_offset = ifd_offsets[ifd_index] as usize;
393            debug_assert_eq!(bytes.len(), ifd_offset);
394
395            let ifd_size = 2 + ifd.entries.len() * 12 + 4;
396            let mut next_data_offset = ifd_offset + ifd_size;
397            let image_offset = next_data_offset as u32;
398            next_data_offset += ifd.image_data.len();
399
400            bytes.extend_from_slice(&le_u16(ifd.entries.len() as u16));
401            let mut deferred = Vec::new();
402            for (tag, ty, count, value) in &ifd.entries {
403                bytes.extend_from_slice(&le_u16(*tag));
404                bytes.extend_from_slice(&le_u16(*ty));
405                bytes.extend_from_slice(&le_u32(*count));
406                if *tag == 273 {
407                    bytes.extend_from_slice(&le_u32(image_offset));
408                } else if value.len() <= 4 {
409                    let mut inline = [0u8; 4];
410                    inline[..value.len()].copy_from_slice(value);
411                    bytes.extend_from_slice(&inline);
412                } else {
413                    bytes.extend_from_slice(&le_u32(next_data_offset as u32));
414                    next_data_offset += value.len();
415                    deferred.push(value.clone());
416                }
417            }
418
419            let next_ifd_offset = ifd_offsets.get(ifd_index + 1).copied().unwrap_or(0);
420            bytes.extend_from_slice(&le_u32(next_ifd_offset));
421            bytes.extend_from_slice(&ifd.image_data);
422            for value in deferred {
423                bytes.extend_from_slice(&value);
424            }
425            debug_assert_eq!(bytes.len(), next_data_offset);
426        }
427
428        bytes
429    }
430
431    fn build_simple_geotiff(pixel_is_point: bool) -> Vec<u8> {
432        let image_data = vec![10u8, 20, 30, 40];
433        let tiepoints = [0.0, 0.0, 0.0, 100.0, 200.0, 0.0];
434        let scales = [2.0, 2.0, 0.0];
435        let geo_keys = if pixel_is_point {
436            vec![
437                1, 1, 0, 3, // header
438                1024, 0, 1, 2, // model type = Geographic
439                1025, 0, 1, 2, // raster type = PixelIsPoint
440                2048, 0, 1, 4326, // EPSG:4326
441            ]
442        } else {
443            vec![
444                1, 1, 0, 2, // header
445                1024, 0, 1, 2, // model type = Geographic
446                2048, 0, 1, 4326, // EPSG:4326
447            ]
448        };
449        let nodata = b"-9999\0".to_vec();
450
451        build_classic_tiff(&[TestIfdSpec {
452            image_data,
453            entries: vec![
454                (256u16, 4u16, 1u32, le_u32(2).to_vec()),
455                (257u16, 4u16, 1u32, le_u32(2).to_vec()),
456                (258u16, 3u16, 1u32, [8, 0, 0, 0].to_vec()),
457                (259u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
458                (273u16, 4u16, 1u32, vec![]),
459                (277u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
460                (278u16, 4u16, 1u32, le_u32(2).to_vec()),
461                (279u16, 4u16, 1u32, le_u32(4).to_vec()),
462                (
463                    33550u16,
464                    12u16,
465                    3u32,
466                    scales.iter().flat_map(|value| le_f64(*value)).collect(),
467                ),
468                (
469                    33922u16,
470                    12u16,
471                    6u32,
472                    tiepoints.iter().flat_map(|value| le_f64(*value)).collect(),
473                ),
474                (
475                    34735u16,
476                    3u16,
477                    geo_keys.len() as u32,
478                    geo_keys.iter().flat_map(|value| le_u16(*value)).collect(),
479                ),
480                (42113u16, 2u16, nodata.len() as u32, nodata),
481            ],
482        }])
483    }
484
485    fn build_geotiff_with_overview() -> Vec<u8> {
486        let base = TestIfdSpec {
487            image_data: vec![10u8, 20, 30, 40],
488            entries: vec![
489                (256u16, 4u16, 1u32, le_u32(2).to_vec()),
490                (257u16, 4u16, 1u32, le_u32(2).to_vec()),
491                (258u16, 3u16, 1u32, [8, 0, 0, 0].to_vec()),
492                (259u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
493                (273u16, 4u16, 1u32, vec![]),
494                (277u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
495                (278u16, 4u16, 1u32, le_u32(2).to_vec()),
496                (279u16, 4u16, 1u32, le_u32(4).to_vec()),
497                (
498                    33550u16,
499                    12u16,
500                    3u32,
501                    [2.0, 2.0, 0.0]
502                        .iter()
503                        .flat_map(|value| le_f64(*value))
504                        .collect(),
505                ),
506                (
507                    33922u16,
508                    12u16,
509                    6u32,
510                    [0.0, 0.0, 0.0, 100.0, 200.0, 0.0]
511                        .iter()
512                        .flat_map(|value| le_f64(*value))
513                        .collect(),
514                ),
515                (
516                    34735u16,
517                    3u16,
518                    12u32,
519                    [1u16, 1, 0, 2, 1024, 0, 1, 2, 2048, 0, 1, 4326]
520                        .iter()
521                        .flat_map(|value| le_u16(*value))
522                        .collect(),
523                ),
524            ],
525        };
526        let overview = TestIfdSpec {
527            image_data: vec![99u8],
528            entries: vec![
529                (254u16, 4u16, 1u32, le_u32(1).to_vec()),
530                (256u16, 4u16, 1u32, le_u32(1).to_vec()),
531                (257u16, 4u16, 1u32, le_u32(1).to_vec()),
532                (258u16, 3u16, 1u32, [8, 0, 0, 0].to_vec()),
533                (259u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
534                (273u16, 4u16, 1u32, vec![]),
535                (277u16, 3u16, 1u32, [1, 0, 0, 0].to_vec()),
536                (278u16, 4u16, 1u32, le_u32(1).to_vec()),
537                (279u16, 4u16, 1u32, le_u32(1).to_vec()),
538            ],
539        };
540
541        build_classic_tiff(&[base, overview])
542    }
543
544    #[test]
545    fn parses_geotiff_metadata_and_reads_raster() {
546        let file = GeoTiffFile::from_bytes(build_simple_geotiff(false)).unwrap();
547        assert_eq!(file.epsg(), Some(4326));
548        assert_eq!(file.width(), 2);
549        assert_eq!(file.height(), 2);
550        assert_eq!(file.band_count(), 1);
551        assert_eq!(file.nodata(), Some("-9999"));
552        assert_eq!(file.geo_bounds(), Some([100.0, 196.0, 104.0, 200.0]));
553
554        let raster = file.read_raster::<u8>().unwrap();
555        assert_eq!(raster.shape(), &[2, 2]);
556        let (values, offset) = raster.into_raw_vec_and_offset();
557        assert_eq!(offset, Some(0));
558        assert_eq!(values, vec![10, 20, 30, 40]);
559    }
560
561    #[test]
562    fn pixel_is_point_metadata_shifts_bounds_to_outer_edges() {
563        let file = GeoTiffFile::from_bytes(build_simple_geotiff(true)).unwrap();
564        assert_eq!(file.geo_bounds(), Some([99.0, 197.0, 103.0, 201.0]));
565
566        let transform = file.transform().unwrap();
567        let (center_x, center_y) = transform.pixel_to_geo(0.5, 0.5);
568        assert_eq!((center_x, center_y), (100.0, 200.0));
569    }
570
571    #[test]
572    fn discovers_reduced_resolution_overviews() {
573        let file = GeoTiffFile::from_bytes(build_geotiff_with_overview()).unwrap();
574        assert_eq!(file.overview_count(), 1);
575        assert_eq!(file.overview_ifd_index(0).unwrap(), 1);
576
577        let overview = file.read_overview::<u8>(0).unwrap();
578        assert_eq!(overview.shape(), &[1, 1]);
579        let (values, offset) = overview.into_raw_vec_and_offset();
580        assert_eq!(offset, Some(0));
581        assert_eq!(values, vec![99]);
582    }
583
584    #[test]
585    fn reads_base_raster_window() {
586        let file = GeoTiffFile::from_bytes(build_simple_geotiff(false)).unwrap();
587        let window = file.read_window::<u8>(1, 0, 1, 2).unwrap();
588        assert_eq!(window.shape(), &[1, 2]);
589        let (values, offset) = window.into_raw_vec_and_offset();
590        assert_eq!(offset, Some(0));
591        assert_eq!(values, vec![30, 40]);
592    }
593
594    #[test]
595    fn reads_overview_window() {
596        let file = GeoTiffFile::from_bytes(build_geotiff_with_overview()).unwrap();
597        let window = file.read_overview_window::<u8>(0, 0, 0, 1, 1).unwrap();
598        assert_eq!(window.shape(), &[1, 1]);
599        let (values, offset) = window.into_raw_vec_and_offset();
600        assert_eq!(offset, Some(0));
601        assert_eq!(values, vec![99]);
602    }
603}
geotiff_reader/lib.rs

geotiff_reader/
lib.rs
