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