npyz

src/data.rs

@@ -1,9 +1,11 @@
 use ndarray::{
-    array, s, Array, Array1, Array2, Array3, ArrayBase, ArrayView2, Ix2, OwnedRepr, ViewRepr,
+    arr2, array, s, Array, Array1, Array2, Array3, ArrayBase, ArrayView2, Axis, Ix2, OwnedRepr,
+    RawDataClone, Slice, ViewRepr,
 };
+use npyz::{npz::NpzArchive, DType, Deserialize, TypeChar};
 use num_traits::{AsPrimitive, FromPrimitive, Num, ToPrimitive};
 use quadtree_rs::{area::AreaBuilder, Quadtree};
-use std::{self, borrow::Borrow, f64::consts::PI};
+use std::{self, borrow::Borrow, f64::consts::PI, io::BufReader, path::Path};
 use thiserror::Error;
 
 type Lon = f64;
@@ -53,9 +55,7 @@ where
     pub data: Array3<T>,
 }
 
-trait MultiDimensionData<T: Num> {
+trait MultiDimensionData {}
-    fn levels(&self, levels: impl IntoIterator<Item = T>) -> Self;
-}
 
 impl<T, Raw> RadarData2d<T, Raw>
 where
@@ -76,7 +76,7 @@ where
         width_rate: f64,
         height_rate: f64,
         filter_len: f64,
-    ) -> Result<Array2<T>, DataError> {
+    ) -> Result<RadarData2d<T, OwnedRepr<T>>, DataError> {
         let width_rate = width_rate.min(1.0);
         let height_rate = height_rate.min(1.0);
         match self.coord_type {
@@ -84,12 +84,31 @@ where
             CoorType::LatLon => {
                 let width_filtered: ArrayBase<ndarray::OwnedRepr<T>, Ix2> =
                     Self::_resample(&self.data, width_rate, filter_len, false);
-                Ok(Self::_resample(
+
-                    &width_filtered,
+                let new_dim1 = Self::_resample(
-                    height_rate,
+                    &Array2::from_shape_vec((1, self.dim1.len()), self.dim1.to_vec()).unwrap(),
+                    width_rate,
                     filter_len,
-                    true,
+                    false,
-                ))
+                )
+                .slice(s![0, ..])
+                .to_owned();
+
+                let new_dim2 = Self::_resample(
+                    &Array2::from_shape_vec((1, self.dim2.len()), self.dim2.to_vec()).unwrap(),
+                    width_rate,
+                    filter_len,
+                    false,
+                )
+                .slice(s![0, ..])
+                .to_owned();
+
+                Ok(RadarData2d {
+                    dim1: new_dim1,
+                    dim2: new_dim2,
+                    data: Self::_resample(&width_filtered, height_rate, filter_len, true),
+                    coord_type: self.coord_type.to_owned(),
+                })
             }
         }
     }
@@ -131,17 +150,20 @@ where
         ];
     }
 
-    fn _resample<'a, R: ndarray::Data<Elem = T>>(
+    fn _resample<'a, V, R: ndarray::Data<Elem = V>>(
         data: &'a ArrayBase<R, Ix2>,
         rate: f64,
         filter_len: f64,
         reverse: bool,
-    ) -> Array2<T> {
+    ) -> Array2<V>
+    where
+        V: Num + Clone + AsPrimitive<f64> + FromPrimitive,
+    {
         let ori_width = if reverse { data.nrows() } else { data.ncols() };
         let ori_height = if reverse { data.ncols() } else { data.nrows() };
 
         let new_width = (ori_width as f64 * rate).ceil() as usize;
-        let mut result: Array2<T> = Array2::zeros((ori_height, new_width));
+        let mut result: Array2<V> = Array2::zeros((ori_height, new_width));
         (0..ori_height).into_iter().for_each(|height| {
             for width in 0..new_width {
                 let center_x = (width as f64 + 0.5) / new_width as f64 * ori_width as f64;
@@ -161,7 +183,7 @@ where
                     filter_sum += weight;
                 }
 
-                result[[height, width]] = T::from_f64(value_sum / filter_sum).unwrap();
+                result[[height, width]] = V::from_f64(value_sum / filter_sum).unwrap();
             }
         });
         result
@@ -180,14 +202,9 @@ fn windowed_sinc(x: f64, y: f64) -> f64 {
     sinc * window
 }
 
-pub struct LevelData<'a, T>
+pub struct LevelData<T: Num + Clone>(Quadtree<i64, RadarData2d<T, OwnedRepr<T>>>);
-where
-    T: Num + Clone,
-{
-    tree: Quadtree<i64, RadarData2d<T, ViewRepr<&'a T>>>,
-}
 
-impl<'a, T> LevelData<'a, T>
+impl<T> LevelData<T>
 where
     T: Num + Clone + AsPrimitive<f64> + FromPrimitive,
 {
@@ -195,7 +212,6 @@ where
         level_data: Vec<RadarData2d<T, OwnedRepr<T>>>,
         level_num: usize,
     ) -> Vec<RadarData2d<T, OwnedRepr<T>>> {
-
         if level_num == 0 {
             return level_data;
         }
@@ -209,13 +225,115 @@ where
         return Self::value(result, level_num - 1);
     }
 
+    fn new(data: RadarData2d<T, OwnedRepr<T>>, level: usize) -> Self {
+        let rate = 1.0 / level as f64;
+        let resampled = data.resample(rate, rate, 2.0).unwrap();
+        let blocks = Self::value(vec![resampled], level);
+        let mut tree: Quadtree<i64, RadarData2d<T, OwnedRepr<T>>> =
+            quadtree_rs::Quadtree::new(blocks.len());
+
+        blocks.into_iter().for_each(|block| {
+            tree.insert(
+                AreaBuilder::default()
+                    .anchor(quadtree_rs::point::Point {
+                        x: *block.dim1.first().unwrap() as i64,
+                        y: *block.dim2.first().unwrap() as i64,
+                    })
+                    .dimensions((block.dim1.len() as i64, block.dim2.len() as i64))
+                    .build()
+                    .unwrap(),
+                block,
+            );
+        });
+
+        Self(tree)
+    }
 }
 
-// impl<T: Num + AsPrimitive<f64> + FromPrimitive> MultiDimensionData<T> for RadarData2d<T> {
+pub fn levels<T>(data: RadarData2d<T, OwnedRepr<T>>, levels: usize) -> Vec<LevelData<T>>
-//     fn levels(&self, levels: usize) -> Self {
+where
-//         let a: Vec<_> = Array1::linspace(0.0, 1.0, levels)
+    T: Num + Clone + AsPrimitive<f64> + FromPrimitive,
-//             .into_iter()
+{
-//             .map(|scale| self.resample(scale, scale, 3.0))
+    (0..levels)
-//             .collect();
+        .into_iter()
-//     }
+        .map(|level| LevelData::new(data.clone(), level))
-// }
+        .collect()
+}
+
+impl<T, Raw> MultiDimensionData for RadarData2d<T, Raw>
+where
+    T: Num + Clone,
+    Raw: ndarray::Data<Elem = T> + Clone + RawDataClone,
+{
+}
+
+pub trait DataLoader<T: MultiDimensionData> {
+    fn load<P: AsRef<Path>>(&self, path: P) -> Result<T, DataError>;
+}
+
+pub struct Npz;
+
+impl Npz {
+    #[inline]
+    fn load_1d<T: Num + Deserialize>(
+        &self,
+        // path: &Path,
+        data: &mut NpzArchive<BufReader<std::fs::File>>,
+        name: &str,
+    ) -> Result<Array1<T>, DataError> {
+        // let mut data = npyz::npz::NpzArchive::open(path)?;
+        let a = data.by_name(name)?.unwrap();
+        let b: Vec<T> = a.into_vec().unwrap();
+
+        Ok(Array1::from_shape_vec(b.len(), b).unwrap())
+    }
+    #[inline]
+    fn load_2d<T: Num + Deserialize>(
+        &self,
+        // path: &Path,
+        data: &mut NpzArchive<BufReader<std::fs::File>>,
+        name: &str,
+    ) -> Result<Array2<T>, DataError> {
+        // let mut data = npyz::npz::NpzArchive::open(path)?;
+        let a = data.by_name(name)?.unwrap();
+
+        let shape = a.shape().to_vec();
+        let b: Vec<T> = a.into_vec().unwrap();
+
+        Ok(Array2::from_shape_vec((shape[0] as usize, shape[1] as usize), b).unwrap())
+    }
+    fn load_3d(&self, data: &mut NpzArchive<BufReader<std::fs::File>>) {}
+}
+
+impl<T> DataLoader<RadarData2d<T, OwnedRepr<T>>> for Npz
+where
+    T: Num + Clone + Deserialize,
+{
+    fn load<P: AsRef<Path>>(&self, path: P) -> Result<RadarData2d<T, OwnedRepr<T>>, DataError> {
+        let mut data: NpzArchive<BufReader<std::fs::File>> = npyz::npz::NpzArchive::open(path)?;
+        if let DType::Plain(type_str) = data.by_name("val")?.unwrap().dtype() {
+            let bytes_num = type_str.num_bytes().unwrap();
+            match type_str.type_char() {
+                TypeChar::Float => match bytes_num {
+                    4 => self.load_2d::<f32>(&mut data, "val")?,
+                    8 => self.load_2d::<f64>(&mut data, "val")?,
+                    _ => {}
+                },
+
+                TypeChar::Int => {}
+
+                _ => {}
+            }
+        }
+        let dim1 = self.load_1d::<f64>(&mut data, "dim1")?;
+        let dim2 = self.load_1d::<f64>(&mut data, "dim2")?;
+        let value = self.load_2d::<T>(&mut data, "val")?;
+
+        Ok(RadarData2d {
+            dim1: dim1,
+            dim2: dim2,
+            data: value,
+            coord_type: CoorType::LatLon,
+        })
+    }
+}