1#![allow(unsafe_code)]
4
5use core::fmt;
75
76#[cfg(not(feature = "std"))]
77use crate::compat::*;
78use crate::error::{OxiGdalError, Result};
79use crate::types::{NoDataValue, RasterDataType};
80
81mod band_iterator;
82pub use band_iterator::{BandIterator, BandRef, MultiBandBuffer};
83
84mod raster_window;
85pub use raster_window::RasterWindow;
86
87pub mod mask;
88pub use mask::Mask;
89
90#[cfg(feature = "std")]
91pub use crate::simd_buffer::{ArenaTile, TileIteratorArena};
92
93#[cfg(feature = "arrow")]
94pub mod arrow_convert;
95
96#[derive(Clone)]
98pub struct RasterBuffer {
99 data: Vec<u8>,
101 width: u64,
103 height: u64,
105 data_type: RasterDataType,
107 nodata: NoDataValue,
109}
110
111impl RasterBuffer {
112 pub fn new(
117 data: Vec<u8>,
118 width: u64,
119 height: u64,
120 data_type: RasterDataType,
121 nodata: NoDataValue,
122 ) -> Result<Self> {
123 let expected_size = width * height * data_type.size_bytes() as u64;
124 if data.len() as u64 != expected_size {
125 return Err(OxiGdalError::InvalidParameter {
126 parameter: "data",
127 message: format!(
128 "Data size mismatch: expected {} bytes for {}x{} {:?}, got {}",
129 expected_size,
130 width,
131 height,
132 data_type,
133 data.len()
134 ),
135 });
136 }
137
138 Ok(Self {
139 data,
140 width,
141 height,
142 data_type,
143 nodata,
144 })
145 }
146
147 #[must_use]
149 pub fn zeros(width: u64, height: u64, data_type: RasterDataType) -> Self {
150 let size = (width * height * data_type.size_bytes() as u64) as usize;
151 Self {
152 data: vec![0u8; size],
153 width,
154 height,
155 data_type,
156 nodata: NoDataValue::None,
157 }
158 }
159
160 #[must_use]
162 pub fn nodata_filled(
163 width: u64,
164 height: u64,
165 data_type: RasterDataType,
166 nodata: NoDataValue,
167 ) -> Self {
168 let mut buffer = Self::zeros(width, height, data_type);
169 buffer.nodata = nodata;
170
171 if let Some(value) = nodata.as_f64() {
173 buffer.fill_value(value);
174 }
175
176 buffer
177 }
178
179 pub fn fill_value(&mut self, value: f64) {
181 match self.data_type {
182 RasterDataType::UInt8 => {
183 let v = value as u8;
184 self.data.fill(v);
185 }
186 RasterDataType::Int8 => {
187 let v = value as i8;
188 self.data.fill(v as u8);
189 }
190 RasterDataType::UInt16 => {
191 let v = (value as u16).to_ne_bytes();
192 for chunk in self.data.chunks_exact_mut(2) {
193 chunk.copy_from_slice(&v);
194 }
195 }
196 RasterDataType::Int16 => {
197 let v = (value as i16).to_ne_bytes();
198 for chunk in self.data.chunks_exact_mut(2) {
199 chunk.copy_from_slice(&v);
200 }
201 }
202 RasterDataType::UInt32 => {
203 let v = (value as u32).to_ne_bytes();
204 for chunk in self.data.chunks_exact_mut(4) {
205 chunk.copy_from_slice(&v);
206 }
207 }
208 RasterDataType::Int32 => {
209 let v = (value as i32).to_ne_bytes();
210 for chunk in self.data.chunks_exact_mut(4) {
211 chunk.copy_from_slice(&v);
212 }
213 }
214 RasterDataType::Float32 => {
215 let v = (value as f32).to_ne_bytes();
216 for chunk in self.data.chunks_exact_mut(4) {
217 chunk.copy_from_slice(&v);
218 }
219 }
220 RasterDataType::Float64 => {
221 let v = value.to_ne_bytes();
222 for chunk in self.data.chunks_exact_mut(8) {
223 chunk.copy_from_slice(&v);
224 }
225 }
226 RasterDataType::UInt64 => {
227 let v = (value as u64).to_ne_bytes();
228 for chunk in self.data.chunks_exact_mut(8) {
229 chunk.copy_from_slice(&v);
230 }
231 }
232 RasterDataType::Int64 => {
233 let v = (value as i64).to_ne_bytes();
234 for chunk in self.data.chunks_exact_mut(8) {
235 chunk.copy_from_slice(&v);
236 }
237 }
238 RasterDataType::CFloat32 => {
239 let real_bytes = (value as f32).to_ne_bytes();
242 let imag_bytes = 0f32.to_ne_bytes();
243 for chunk in self.data.chunks_exact_mut(8) {
244 chunk[..4].copy_from_slice(&real_bytes);
245 chunk[4..].copy_from_slice(&imag_bytes);
246 }
247 }
248 RasterDataType::CFloat64 => {
249 let real_bytes = value.to_ne_bytes();
252 let imag_bytes = 0f64.to_ne_bytes();
253 for chunk in self.data.chunks_exact_mut(16) {
254 chunk[..8].copy_from_slice(&real_bytes);
255 chunk[8..].copy_from_slice(&imag_bytes);
256 }
257 }
258 }
259 }
260
261 #[must_use]
263 pub const fn width(&self) -> u64 {
264 self.width
265 }
266
267 #[must_use]
269 pub const fn height(&self) -> u64 {
270 self.height
271 }
272
273 #[must_use]
275 pub const fn data_type(&self) -> RasterDataType {
276 self.data_type
277 }
278
279 #[must_use]
281 pub const fn nodata(&self) -> NoDataValue {
282 self.nodata
283 }
284
285 #[must_use]
287 pub const fn pixel_count(&self) -> u64 {
288 self.width * self.height
289 }
290
291 #[must_use]
293 pub fn as_bytes(&self) -> &[u8] {
294 &self.data
295 }
296
297 pub fn as_bytes_mut(&mut self) -> &mut [u8] {
299 &mut self.data
300 }
301
302 #[must_use]
304 pub fn into_bytes(self) -> Vec<u8> {
305 self.data
306 }
307
308 pub fn from_typed_vec<T: Copy + 'static>(
319 width: usize,
320 height: usize,
321 data: Vec<T>,
322 data_type: RasterDataType,
323 ) -> Result<Self> {
324 let expected_pixels = width * height;
325 if data.len() != expected_pixels {
326 return Err(OxiGdalError::InvalidParameter {
327 parameter: "data",
328 message: format!(
329 "Data length mismatch: expected {} pixels for {}x{}, got {}",
330 expected_pixels,
331 width,
332 height,
333 data.len()
334 ),
335 });
336 }
337
338 let type_size = core::mem::size_of::<T>();
340 let expected_type_size = data_type.size_bytes();
341 if type_size != expected_type_size {
342 return Err(OxiGdalError::InvalidParameter {
343 parameter: "data_type",
344 message: format!(
345 "Type size mismatch: provided type has {} bytes, {:?} expects {} bytes",
346 type_size, data_type, expected_type_size
347 ),
348 });
349 }
350
351 let byte_data: Vec<u8> = data
352 .iter()
353 .flat_map(|v| {
354 let ptr = v as *const T as *const u8;
356 unsafe { core::slice::from_raw_parts(ptr, type_size) }.to_vec()
357 })
358 .collect();
359
360 Self::new(
361 byte_data,
362 width as u64,
363 height as u64,
364 data_type,
365 NoDataValue::None,
366 )
367 }
368
369 pub fn as_slice<T: Copy + 'static>(&self) -> Result<&[T]> {
377 let type_size = core::mem::size_of::<T>();
378 let expected_size = self.data_type.size_bytes();
379
380 if type_size != expected_size {
381 return Err(OxiGdalError::InvalidParameter {
382 parameter: "T",
383 message: format!(
384 "Type size mismatch: requested type has {} bytes, buffer contains {:?} ({} bytes)",
385 type_size, self.data_type, expected_size
386 ),
387 });
388 }
389
390 let pixel_count = (self.width * self.height) as usize;
391 let slice =
394 unsafe { core::slice::from_raw_parts(self.data.as_ptr() as *const T, pixel_count) };
395 Ok(slice)
396 }
397
398 pub fn as_slice_mut<T: Copy + 'static>(&mut self) -> Result<&mut [T]> {
406 let type_size = core::mem::size_of::<T>();
407 let expected_size = self.data_type.size_bytes();
408
409 if type_size != expected_size {
410 return Err(OxiGdalError::InvalidParameter {
411 parameter: "T",
412 message: format!(
413 "Type size mismatch: requested type has {} bytes, buffer contains {:?} ({} bytes)",
414 type_size, self.data_type, expected_size
415 ),
416 });
417 }
418
419 let pixel_count = (self.width * self.height) as usize;
420 let slice = unsafe {
423 core::slice::from_raw_parts_mut(self.data.as_mut_ptr() as *mut T, pixel_count)
424 };
425 Ok(slice)
426 }
427
428 pub fn get_pixel(&self, x: u64, y: u64) -> Result<f64> {
433 if x >= self.width || y >= self.height {
434 return Err(OxiGdalError::OutOfBounds {
435 message: format!(
436 "Pixel ({}, {}) out of bounds for {}x{} buffer",
437 x, y, self.width, self.height
438 ),
439 });
440 }
441
442 let pixel_size = self.data_type.size_bytes();
443 let offset = (y * self.width + x) as usize * pixel_size;
444
445 let value = match self.data_type {
446 RasterDataType::UInt8 => f64::from(self.data[offset]),
447 RasterDataType::Int8 => f64::from(self.data[offset] as i8),
448 RasterDataType::UInt16 => {
449 let bytes: [u8; 2] = self.data[offset..offset + 2].try_into().map_err(|_| {
450 OxiGdalError::Internal {
451 message: "Invalid slice length".to_string(),
452 }
453 })?;
454 f64::from(u16::from_ne_bytes(bytes))
455 }
456 RasterDataType::Int16 => {
457 let bytes: [u8; 2] = self.data[offset..offset + 2].try_into().map_err(|_| {
458 OxiGdalError::Internal {
459 message: "Invalid slice length".to_string(),
460 }
461 })?;
462 f64::from(i16::from_ne_bytes(bytes))
463 }
464 RasterDataType::UInt32 => {
465 let bytes: [u8; 4] = self.data[offset..offset + 4].try_into().map_err(|_| {
466 OxiGdalError::Internal {
467 message: "Invalid slice length".to_string(),
468 }
469 })?;
470 f64::from(u32::from_ne_bytes(bytes))
471 }
472 RasterDataType::Int32 => {
473 let bytes: [u8; 4] = self.data[offset..offset + 4].try_into().map_err(|_| {
474 OxiGdalError::Internal {
475 message: "Invalid slice length".to_string(),
476 }
477 })?;
478 f64::from(i32::from_ne_bytes(bytes))
479 }
480 RasterDataType::Float32 => {
481 let bytes: [u8; 4] = self.data[offset..offset + 4].try_into().map_err(|_| {
482 OxiGdalError::Internal {
483 message: "Invalid slice length".to_string(),
484 }
485 })?;
486 f64::from(f32::from_ne_bytes(bytes))
487 }
488 RasterDataType::Float64 => {
489 let bytes: [u8; 8] = self.data[offset..offset + 8].try_into().map_err(|_| {
490 OxiGdalError::Internal {
491 message: "Invalid slice length".to_string(),
492 }
493 })?;
494 f64::from_ne_bytes(bytes)
495 }
496 RasterDataType::UInt64 => {
497 let bytes: [u8; 8] = self.data[offset..offset + 8].try_into().map_err(|_| {
498 OxiGdalError::Internal {
499 message: "Invalid slice length".to_string(),
500 }
501 })?;
502 u64::from_ne_bytes(bytes) as f64
503 }
504 RasterDataType::Int64 => {
505 let bytes: [u8; 8] = self.data[offset..offset + 8].try_into().map_err(|_| {
506 OxiGdalError::Internal {
507 message: "Invalid slice length".to_string(),
508 }
509 })?;
510 i64::from_ne_bytes(bytes) as f64
511 }
512 RasterDataType::CFloat32 => {
513 let bytes: [u8; 4] = self.data[offset..offset + 4].try_into().map_err(|_| {
515 OxiGdalError::Internal {
516 message: "Invalid slice length".to_string(),
517 }
518 })?;
519 f64::from(f32::from_ne_bytes(bytes))
520 }
521 RasterDataType::CFloat64 => {
522 let bytes: [u8; 8] = self.data[offset..offset + 8].try_into().map_err(|_| {
524 OxiGdalError::Internal {
525 message: "Invalid slice length".to_string(),
526 }
527 })?;
528 f64::from_ne_bytes(bytes)
529 }
530 };
531
532 Ok(value)
533 }
534
535 pub fn set_pixel(&mut self, x: u64, y: u64, value: f64) -> Result<()> {
540 if x >= self.width || y >= self.height {
541 return Err(OxiGdalError::OutOfBounds {
542 message: format!(
543 "Pixel ({}, {}) out of bounds for {}x{} buffer",
544 x, y, self.width, self.height
545 ),
546 });
547 }
548
549 let pixel_size = self.data_type.size_bytes();
550 let offset = (y * self.width + x) as usize * pixel_size;
551
552 match self.data_type {
553 RasterDataType::UInt8 => {
554 self.data[offset] = value as u8;
555 }
556 RasterDataType::Int8 => {
557 self.data[offset] = (value as i8) as u8;
558 }
559 RasterDataType::UInt16 => {
560 let bytes = (value as u16).to_ne_bytes();
561 self.data[offset..offset + 2].copy_from_slice(&bytes);
562 }
563 RasterDataType::Int16 => {
564 let bytes = (value as i16).to_ne_bytes();
565 self.data[offset..offset + 2].copy_from_slice(&bytes);
566 }
567 RasterDataType::UInt32 => {
568 let bytes = (value as u32).to_ne_bytes();
569 self.data[offset..offset + 4].copy_from_slice(&bytes);
570 }
571 RasterDataType::Int32 => {
572 let bytes = (value as i32).to_ne_bytes();
573 self.data[offset..offset + 4].copy_from_slice(&bytes);
574 }
575 RasterDataType::Float32 => {
576 let bytes = (value as f32).to_ne_bytes();
577 self.data[offset..offset + 4].copy_from_slice(&bytes);
578 }
579 RasterDataType::Float64 => {
580 let bytes = value.to_ne_bytes();
581 self.data[offset..offset + 8].copy_from_slice(&bytes);
582 }
583 RasterDataType::UInt64 => {
584 let bytes = (value as u64).to_ne_bytes();
585 self.data[offset..offset + 8].copy_from_slice(&bytes);
586 }
587 RasterDataType::Int64 => {
588 let bytes = (value as i64).to_ne_bytes();
589 self.data[offset..offset + 8].copy_from_slice(&bytes);
590 }
591 RasterDataType::CFloat32 => {
592 let bytes = (value as f32).to_ne_bytes();
594 self.data[offset..offset + 4].copy_from_slice(&bytes);
595 }
596 RasterDataType::CFloat64 => {
597 let bytes = value.to_ne_bytes();
599 self.data[offset..offset + 8].copy_from_slice(&bytes);
600 }
601 }
602
603 Ok(())
604 }
605
606 pub fn get_u8(&self, x: u64, y: u64) -> Result<u8> {
613 self.check_bounds(x, y)?;
614 self.check_type(RasterDataType::UInt8)?;
615 let offset = (y * self.width + x) as usize;
616 Ok(self.data[offset])
617 }
618
619 pub fn get_i8(&self, x: u64, y: u64) -> Result<i8> {
624 self.check_bounds(x, y)?;
625 self.check_type(RasterDataType::Int8)?;
626 let offset = (y * self.width + x) as usize;
627 Ok(self.data[offset] as i8)
628 }
629
630 pub fn get_u16(&self, x: u64, y: u64) -> Result<u16> {
635 self.check_bounds(x, y)?;
636 self.check_type(RasterDataType::UInt16)?;
637 let offset = (y * self.width + x) as usize * 2;
638 let bytes: [u8; 2] =
639 self.data[offset..offset + 2]
640 .try_into()
641 .map_err(|_| OxiGdalError::Internal {
642 message: "Invalid slice length".to_string(),
643 })?;
644 Ok(u16::from_ne_bytes(bytes))
645 }
646
647 pub fn get_i16(&self, x: u64, y: u64) -> Result<i16> {
652 self.check_bounds(x, y)?;
653 self.check_type(RasterDataType::Int16)?;
654 let offset = (y * self.width + x) as usize * 2;
655 let bytes: [u8; 2] =
656 self.data[offset..offset + 2]
657 .try_into()
658 .map_err(|_| OxiGdalError::Internal {
659 message: "Invalid slice length".to_string(),
660 })?;
661 Ok(i16::from_ne_bytes(bytes))
662 }
663
664 pub fn get_u32(&self, x: u64, y: u64) -> Result<u32> {
669 self.check_bounds(x, y)?;
670 self.check_type(RasterDataType::UInt32)?;
671 let offset = (y * self.width + x) as usize * 4;
672 let bytes: [u8; 4] =
673 self.data[offset..offset + 4]
674 .try_into()
675 .map_err(|_| OxiGdalError::Internal {
676 message: "Invalid slice length".to_string(),
677 })?;
678 Ok(u32::from_ne_bytes(bytes))
679 }
680
681 pub fn get_i32(&self, x: u64, y: u64) -> Result<i32> {
686 self.check_bounds(x, y)?;
687 self.check_type(RasterDataType::Int32)?;
688 let offset = (y * self.width + x) as usize * 4;
689 let bytes: [u8; 4] =
690 self.data[offset..offset + 4]
691 .try_into()
692 .map_err(|_| OxiGdalError::Internal {
693 message: "Invalid slice length".to_string(),
694 })?;
695 Ok(i32::from_ne_bytes(bytes))
696 }
697
698 pub fn get_u64(&self, x: u64, y: u64) -> Result<u64> {
703 self.check_bounds(x, y)?;
704 self.check_type(RasterDataType::UInt64)?;
705 let offset = (y * self.width + x) as usize * 8;
706 let bytes: [u8; 8] =
707 self.data[offset..offset + 8]
708 .try_into()
709 .map_err(|_| OxiGdalError::Internal {
710 message: "Invalid slice length".to_string(),
711 })?;
712 Ok(u64::from_ne_bytes(bytes))
713 }
714
715 pub fn get_i64(&self, x: u64, y: u64) -> Result<i64> {
720 self.check_bounds(x, y)?;
721 self.check_type(RasterDataType::Int64)?;
722 let offset = (y * self.width + x) as usize * 8;
723 let bytes: [u8; 8] =
724 self.data[offset..offset + 8]
725 .try_into()
726 .map_err(|_| OxiGdalError::Internal {
727 message: "Invalid slice length".to_string(),
728 })?;
729 Ok(i64::from_ne_bytes(bytes))
730 }
731
732 pub fn get_f32(&self, x: u64, y: u64) -> Result<f32> {
737 self.check_bounds(x, y)?;
738 self.check_type(RasterDataType::Float32)?;
739 let offset = (y * self.width + x) as usize * 4;
740 let bytes: [u8; 4] =
741 self.data[offset..offset + 4]
742 .try_into()
743 .map_err(|_| OxiGdalError::Internal {
744 message: "Invalid slice length".to_string(),
745 })?;
746 Ok(f32::from_ne_bytes(bytes))
747 }
748
749 pub fn get_f64(&self, x: u64, y: u64) -> Result<f64> {
754 self.check_bounds(x, y)?;
755 self.check_type(RasterDataType::Float64)?;
756 let offset = (y * self.width + x) as usize * 8;
757 let bytes: [u8; 8] =
758 self.data[offset..offset + 8]
759 .try_into()
760 .map_err(|_| OxiGdalError::Internal {
761 message: "Invalid slice length".to_string(),
762 })?;
763 Ok(f64::from_ne_bytes(bytes))
764 }
765
766 pub fn set_u8(&mut self, x: u64, y: u64, value: u8) -> Result<()> {
773 self.check_bounds(x, y)?;
774 self.check_type(RasterDataType::UInt8)?;
775 let offset = (y * self.width + x) as usize;
776 self.data[offset] = value;
777 Ok(())
778 }
779
780 pub fn set_f32(&mut self, x: u64, y: u64, value: f32) -> Result<()> {
785 self.check_bounds(x, y)?;
786 self.check_type(RasterDataType::Float32)?;
787 let offset = (y * self.width + x) as usize * 4;
788 self.data[offset..offset + 4].copy_from_slice(&value.to_ne_bytes());
789 Ok(())
790 }
791
792 pub fn set_f64(&mut self, x: u64, y: u64, value: f64) -> Result<()> {
797 self.check_bounds(x, y)?;
798 self.check_type(RasterDataType::Float64)?;
799 let offset = (y * self.width + x) as usize * 8;
800 self.data[offset..offset + 8].copy_from_slice(&value.to_ne_bytes());
801 Ok(())
802 }
803
804 pub fn row_slice<T: Copy + 'static>(&self, y: u64) -> Result<&[T]> {
811 if y >= self.height {
812 return Err(OxiGdalError::OutOfBounds {
813 message: format!("Row {} out of bounds for height {}", y, self.height),
814 });
815 }
816 let type_size = core::mem::size_of::<T>();
817 let expected_size = self.data_type.size_bytes();
818 if type_size != expected_size {
819 return Err(OxiGdalError::InvalidParameter {
820 parameter: "T",
821 message: format!(
822 "Type size {} doesn't match {:?} size {}",
823 type_size, self.data_type, expected_size
824 ),
825 });
826 }
827 let row_start = (y * self.width) as usize * expected_size;
828 let row_end = row_start + self.width as usize * expected_size;
829 let slice = unsafe {
831 core::slice::from_raw_parts(
832 self.data[row_start..row_end].as_ptr() as *const T,
833 self.width as usize,
834 )
835 };
836 Ok(slice)
837 }
838
839 pub fn window(&self, x: u64, y: u64, width: u64, height: u64) -> Result<Self> {
844 if x + width > self.width || y + height > self.height {
845 return Err(OxiGdalError::OutOfBounds {
846 message: format!(
847 "Window ({},{}) {}x{} exceeds buffer {}x{}",
848 x, y, width, height, self.width, self.height
849 ),
850 });
851 }
852 let pixel_size = self.data_type.size_bytes();
853 let row_bytes = width as usize * pixel_size;
854 let mut data = Vec::with_capacity(height as usize * row_bytes);
855 for row in y..y + height {
856 let src_start = (row * self.width + x) as usize * pixel_size;
857 data.extend_from_slice(&self.data[src_start..src_start + row_bytes]);
858 }
859 Self::new(data, width, height, self.data_type, self.nodata)
860 }
861
862 fn check_bounds(&self, x: u64, y: u64) -> Result<()> {
865 if x >= self.width || y >= self.height {
866 return Err(OxiGdalError::OutOfBounds {
867 message: format!(
868 "Pixel ({}, {}) out of bounds for {}x{} buffer",
869 x, y, self.width, self.height
870 ),
871 });
872 }
873 Ok(())
874 }
875
876 fn check_type(&self, expected: RasterDataType) -> Result<()> {
877 if self.data_type != expected {
878 return Err(OxiGdalError::InvalidParameter {
879 parameter: "data_type",
880 message: format!(
881 "Buffer contains {:?} data, requested {:?}",
882 self.data_type, expected
883 ),
884 });
885 }
886 Ok(())
887 }
888
889 #[must_use]
891 pub fn is_nodata(&self, value: f64) -> bool {
892 match self.nodata.as_f64() {
893 Some(nd) => {
894 if nd.is_nan() && value.is_nan() {
895 true
896 } else {
897 (nd - value).abs() < f64::EPSILON
898 }
899 }
900 None => false,
901 }
902 }
903
904 pub fn convert_to(&self, target_type: RasterDataType) -> Result<Self> {
909 if target_type == self.data_type {
910 return Ok(self.clone());
911 }
912
913 let mut result = Self::zeros(self.width, self.height, target_type);
914 result.nodata = self.nodata;
915
916 for y in 0..self.height {
917 for x in 0..self.width {
918 let value = self.get_pixel(x, y)?;
919 result.set_pixel(x, y, value)?;
920 }
921 }
922
923 Ok(result)
924 }
925
926 pub fn compute_statistics(&self) -> Result<BufferStatistics> {
928 let mut min = f64::MAX;
929 let mut max = f64::MIN;
930 let mut sum = 0.0;
931 let mut sum_sq = 0.0;
932 let mut valid_count = 0u64;
933
934 for y in 0..self.height {
935 for x in 0..self.width {
936 let value = self.get_pixel(x, y)?;
937 if !self.is_nodata(value) && value.is_finite() {
938 min = min.min(value);
939 max = max.max(value);
940 sum += value;
941 sum_sq += value * value;
942 valid_count += 1;
943 }
944 }
945 }
946
947 if valid_count == 0 {
948 return Ok(BufferStatistics {
949 min: f64::NAN,
950 max: f64::NAN,
951 mean: f64::NAN,
952 std_dev: f64::NAN,
953 valid_count: 0,
954 histogram: None,
955 });
956 }
957
958 let mean = sum / valid_count as f64;
959 let variance = (sum_sq / valid_count as f64) - (mean * mean);
960 let std_dev = variance.sqrt();
961
962 Ok(BufferStatistics {
963 min,
964 max,
965 mean,
966 std_dev,
967 valid_count,
968 histogram: None,
969 })
970 }
971
972 pub fn compute_statistics_with_histogram(&self, bin_count: usize) -> Result<BufferStatistics> {
993 if bin_count == 0 {
994 return Err(OxiGdalError::InvalidParameter {
995 parameter: "bin_count",
996 message: "bin_count must be at least 1".to_string(),
997 });
998 }
999
1000 let mut min = f64::MAX;
1002 let mut max = f64::MIN;
1003 let mut sum = 0.0f64;
1004 let mut sum_sq = 0.0f64;
1005 let mut valid_count = 0u64;
1006 let total_pixels = (self.width * self.height) as usize;
1008 let mut valid_values: Vec<f64> = Vec::with_capacity(total_pixels);
1009
1010 for y in 0..self.height {
1011 for x in 0..self.width {
1012 let value = self.get_pixel(x, y)?;
1013 if !self.is_nodata(value) && value.is_finite() {
1014 min = min.min(value);
1015 max = max.max(value);
1016 sum += value;
1017 sum_sq += value * value;
1018 valid_count += 1;
1019 valid_values.push(value);
1020 }
1021 }
1022 }
1023
1024 let mut bins = vec![0u64; bin_count];
1026
1027 if valid_count > 0 {
1028 let range = max - min;
1029 if range == 0.0 {
1030 bins[0] = valid_count;
1032 } else {
1033 for v in &valid_values {
1034 let bin_idx = (((v - min) / range) * bin_count as f64).floor() as usize;
1035 let bin_idx = bin_idx.min(bin_count - 1);
1037 bins[bin_idx] += 1;
1038 }
1039 }
1040 }
1041
1042 if valid_count == 0 {
1043 return Ok(BufferStatistics {
1044 min: f64::NAN,
1045 max: f64::NAN,
1046 mean: f64::NAN,
1047 std_dev: f64::NAN,
1048 valid_count: 0,
1049 histogram: Some(bins),
1050 });
1051 }
1052
1053 let mean = sum / valid_count as f64;
1054 let variance = (sum_sq / valid_count as f64) - (mean * mean);
1055 let std_dev = variance.max(0.0).sqrt();
1056
1057 Ok(BufferStatistics {
1058 min,
1059 max,
1060 mean,
1061 std_dev,
1062 valid_count,
1063 histogram: Some(bins),
1064 })
1065 }
1066}
1067
1068impl fmt::Debug for RasterBuffer {
1069 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1070 f.debug_struct("RasterBuffer")
1071 .field("width", &self.width)
1072 .field("height", &self.height)
1073 .field("data_type", &self.data_type)
1074 .field("nodata", &self.nodata)
1075 .field("bytes", &self.data.len())
1076 .finish()
1077 }
1078}
1079
1080#[derive(Debug, Clone, PartialEq)]
1085pub struct BufferStatistics {
1086 pub min: f64,
1088 pub max: f64,
1090 pub mean: f64,
1092 pub std_dev: f64,
1094 pub valid_count: u64,
1096 pub histogram: Option<Vec<u64>>,
1101}
1102
1103#[cfg(feature = "arrow")]
1104mod arrow_support {
1105 use arrow_array::{Array, Float64Array};
1108
1109 use super::{OxiGdalError, RasterBuffer, Result};
1110
1111 impl RasterBuffer {
1112 pub fn from_arrow_array<A: Array>(_array: &A, _width: u64, _height: u64) -> Result<Self> {
1117 Err(OxiGdalError::NotSupported {
1120 operation: "Arrow array conversion".to_string(),
1121 })
1122 }
1123
1124 pub fn to_float64_array(&self) -> Result<Float64Array> {
1126 let mut values = Vec::with_capacity(self.pixel_count() as usize);
1127 for y in 0..self.height {
1128 for x in 0..self.width {
1129 values.push(self.get_pixel(x, y)?);
1130 }
1131 }
1132 Ok(Float64Array::from(values))
1133 }
1134 }
1135}
1136
1137#[cfg(test)]
1138mod tests {
1139 #![allow(clippy::expect_used)]
1140
1141 use super::*;
1142
1143 #[test]
1144 fn test_buffer_creation() {
1145 let buffer = RasterBuffer::zeros(100, 100, RasterDataType::UInt8);
1146 assert_eq!(buffer.width(), 100);
1147 assert_eq!(buffer.height(), 100);
1148 assert_eq!(buffer.pixel_count(), 10_000);
1149 assert_eq!(buffer.as_bytes().len(), 10_000);
1150 }
1151
1152 #[test]
1153 fn test_pixel_access() {
1154 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Float32);
1155
1156 buffer.set_pixel(5, 5, 42.0).expect("set should work");
1157 let value = buffer.get_pixel(5, 5).expect("get should work");
1158 assert!((value - 42.0).abs() < f64::EPSILON);
1159
1160 assert!(buffer.get_pixel(100, 0).is_err());
1162 assert!(buffer.set_pixel(0, 100, 0.0).is_err());
1163 }
1164
1165 #[test]
1166 fn test_nodata() {
1167 let buffer = RasterBuffer::nodata_filled(
1168 10,
1169 10,
1170 RasterDataType::Float32,
1171 NoDataValue::Float(-9999.0),
1172 );
1173
1174 assert!(buffer.is_nodata(-9999.0));
1175 assert!(!buffer.is_nodata(0.0));
1176 }
1177
1178 #[test]
1179 fn test_statistics() {
1180 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Float32);
1181
1182 for y in 0..10 {
1184 for x in 0..10 {
1185 let value = (y * 10 + x) as f64;
1186 buffer.set_pixel(x, y, value).expect("set should work");
1187 }
1188 }
1189
1190 let stats = buffer.compute_statistics().expect("stats should work");
1191 assert!((stats.min - 0.0).abs() < f64::EPSILON);
1192 assert!((stats.max - 99.0).abs() < f64::EPSILON);
1193 assert!((stats.mean - 49.5).abs() < 0.01);
1194 assert_eq!(stats.valid_count, 100);
1195 }
1196
1197 #[test]
1198 fn test_data_validation() {
1199 let result = RasterBuffer::new(
1201 vec![0u8; 100],
1202 10,
1203 10,
1204 RasterDataType::UInt16, NoDataValue::None,
1206 );
1207 assert!(result.is_err());
1208
1209 let result = RasterBuffer::new(
1211 vec![0u8; 200],
1212 10,
1213 10,
1214 RasterDataType::UInt16,
1215 NoDataValue::None,
1216 );
1217 assert!(result.is_ok());
1218 }
1219
1220 #[test]
1221 fn test_typed_get_u8() {
1222 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt8);
1223 buffer.set_pixel(3, 4, 200.0).expect("set should work");
1224 assert_eq!(buffer.get_u8(3, 4).expect("get_u8"), 200);
1225 assert!(buffer.get_f32(3, 4).is_err());
1227 }
1228
1229 #[test]
1230 fn test_typed_get_i8() {
1231 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Int8);
1232 buffer.set_pixel(0, 0, -42.0).expect("set should work");
1233 assert_eq!(buffer.get_i8(0, 0).expect("get_i8"), -42);
1234 }
1235
1236 #[test]
1237 fn test_typed_get_u16() {
1238 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt16);
1239 buffer.set_pixel(5, 5, 60000.0).expect("set should work");
1240 assert_eq!(buffer.get_u16(5, 5).expect("get_u16"), 60000);
1241 }
1242
1243 #[test]
1244 fn test_typed_get_i16() {
1245 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Int16);
1246 buffer.set_pixel(2, 3, -1234.0).expect("set should work");
1247 assert_eq!(buffer.get_i16(2, 3).expect("get_i16"), -1234);
1248 }
1249
1250 #[test]
1251 fn test_typed_get_u32() {
1252 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt32);
1253 buffer.set_pixel(0, 0, 100_000.0).expect("set should work");
1254 assert_eq!(buffer.get_u32(0, 0).expect("get_u32"), 100_000);
1255 }
1256
1257 #[test]
1258 fn test_typed_get_i32() {
1259 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Int32);
1260 buffer.set_pixel(1, 1, -50_000.0).expect("set should work");
1261 assert_eq!(buffer.get_i32(1, 1).expect("get_i32"), -50_000);
1262 }
1263
1264 #[test]
1265 fn test_typed_get_u64() {
1266 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt64);
1267 buffer
1268 .set_pixel(0, 0, 1_000_000.0)
1269 .expect("set should work");
1270 assert_eq!(buffer.get_u64(0, 0).expect("get_u64"), 1_000_000);
1271 }
1272
1273 #[test]
1274 fn test_typed_get_i64() {
1275 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Int64);
1276 buffer.set_pixel(0, 0, -999_999.0).expect("set should work");
1277 assert_eq!(buffer.get_i64(0, 0).expect("get_i64"), -999_999);
1278 }
1279
1280 #[test]
1281 fn test_typed_get_f32() {
1282 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Float32);
1283 buffer
1284 .set_f32(7, 8, core::f32::consts::PI)
1285 .expect("set_f32 should work");
1286 let val = buffer.get_f32(7, 8).expect("get_f32");
1287 assert!((val - core::f32::consts::PI).abs() < 1e-5);
1288 }
1289
1290 #[test]
1291 fn test_typed_get_f64() {
1292 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::Float64);
1293 buffer
1294 .set_f64(9, 9, core::f64::consts::E)
1295 .expect("set_f64 should work");
1296 let val = buffer.get_f64(9, 9).expect("get_f64");
1297 assert!((val - core::f64::consts::E).abs() < 1e-9);
1298 }
1299
1300 #[test]
1301 fn test_typed_set_u8() {
1302 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt8);
1303 buffer.set_u8(0, 0, 255).expect("set_u8 should work");
1304 assert_eq!(buffer.get_u8(0, 0).expect("get_u8"), 255);
1305 }
1306
1307 #[test]
1308 fn test_typed_out_of_bounds() {
1309 let buffer = RasterBuffer::zeros(10, 10, RasterDataType::Float32);
1310 assert!(buffer.get_f32(10, 0).is_err());
1311 assert!(buffer.get_f32(0, 10).is_err());
1312 }
1313
1314 #[test]
1315 fn test_typed_wrong_type() {
1316 let buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt8);
1317 assert!(buffer.get_f32(0, 0).is_err());
1318 assert!(buffer.get_u16(0, 0).is_err());
1319 assert!(buffer.get_i32(0, 0).is_err());
1320 assert!(buffer.get_f64(0, 0).is_err());
1321 }
1322
1323 #[test]
1324 fn test_row_slice() {
1325 let mut buffer = RasterBuffer::zeros(5, 3, RasterDataType::Float32);
1326 for x in 0..5 {
1327 buffer
1328 .set_pixel(x, 1, (x + 10) as f64)
1329 .expect("set should work");
1330 }
1331 let row: &[f32] = buffer.row_slice(1).expect("row_slice should work");
1332 assert_eq!(row.len(), 5);
1333 assert!((row[0] - 10.0).abs() < 1e-5);
1334 assert!((row[4] - 14.0).abs() < 1e-5);
1335 }
1336
1337 #[test]
1338 fn test_row_slice_out_of_bounds() {
1339 let buffer = RasterBuffer::zeros(5, 3, RasterDataType::Float32);
1340 assert!(buffer.row_slice::<f32>(3).is_err());
1341 }
1342
1343 #[test]
1344 fn test_window() {
1345 let mut buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt8);
1346 buffer.set_pixel(5, 5, 42.0).expect("set should work");
1347 buffer.set_pixel(6, 6, 99.0).expect("set should work");
1348
1349 let win = buffer.window(4, 4, 4, 4).expect("window should work");
1350 assert_eq!(win.width(), 4);
1351 assert_eq!(win.height(), 4);
1352 let val = win.get_pixel(1, 1).expect("get should work");
1354 assert!((val - 42.0).abs() < f64::EPSILON);
1355 let val = win.get_pixel(2, 2).expect("get should work");
1357 assert!((val - 99.0).abs() < f64::EPSILON);
1358 }
1359
1360 #[test]
1361 fn test_window_out_of_bounds() {
1362 let buffer = RasterBuffer::zeros(10, 10, RasterDataType::UInt8);
1363 assert!(buffer.window(8, 8, 4, 4).is_err());
1364 }
1365
1366 #[test]
1367 fn test_fill_value_cfloat32_writes_real_zero_imag() {
1368 let fill: f64 = 1.23456789;
1370 let mut buf = RasterBuffer::zeros(2, 2, RasterDataType::CFloat32);
1371 buf.fill_value(fill);
1372 let v = buf.get_pixel(0, 0).expect("pixel access");
1374 assert!(
1375 (v - fill).abs() < 1e-5,
1376 "real part should be {fill}, got {v}"
1377 );
1378 let raw = buf.as_bytes();
1380 let imag_bytes: [u8; 4] = raw[4..8].try_into().expect("slice");
1381 let imag = f32::from_ne_bytes(imag_bytes);
1382 assert_eq!(imag, 0.0f32, "imaginary part should be 0.0");
1383 }
1384
1385 #[test]
1386 fn test_fill_value_cfloat64_writes_real_zero_imag() {
1387 let fill: f64 = 9.876_543_21_f64;
1389 let mut buf = RasterBuffer::zeros(2, 2, RasterDataType::CFloat64);
1390 buf.fill_value(fill);
1391 let v = buf.get_pixel(0, 0).expect("pixel access");
1392 assert!((v - fill).abs() < 1e-9, "real part mismatch, got {v}");
1393 let raw = buf.as_bytes();
1395 let imag_bytes: [u8; 8] = raw[8..16].try_into().expect("slice");
1396 let imag = f64::from_ne_bytes(imag_bytes);
1397 assert_eq!(imag, 0.0f64, "imaginary part should be 0.0");
1398 }
1399}