From d6f1341b172ef0695593e430fd016e59f49d211f Mon Sep 17 00:00:00 2001 From: sleptworld Date: Wed, 7 Jun 2023 01:11:43 +0800 Subject: [PATCH] npyz --- src/data.rs | 182 +++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 150 insertions(+), 32 deletions(-) diff --git a/src/data.rs b/src/data.rs index 1a523ad..6e24f42 100644 --- a/src/data.rs +++ b/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, } -trait MultiDimensionData { - fn levels(&self, levels: impl IntoIterator) -> Self; -} +trait MultiDimensionData {} impl RadarData2d where @@ -76,7 +76,7 @@ where width_rate: f64, height_rate: f64, filter_len: f64, - ) -> Result, DataError> { + ) -> Result>, 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, Ix2> = Self::_resample(&self.data, width_rate, filter_len, false); - Ok(Self::_resample( - &width_filtered, - height_rate, + + let new_dim1 = Self::_resample( + &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>( + fn _resample<'a, V, R: ndarray::Data>( data: &'a ArrayBase, rate: f64, filter_len: f64, reverse: bool, - ) -> Array2 { + ) -> Array2 + where + V: Num + Clone + AsPrimitive + 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 = Array2::zeros((ori_height, new_width)); + let mut result: Array2 = 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> -where - T: Num + Clone, -{ - tree: Quadtree>>, -} +pub struct LevelData(Quadtree>>); -impl<'a, T> LevelData<'a, T> +impl LevelData where T: Num + Clone + AsPrimitive + FromPrimitive, { @@ -195,7 +212,6 @@ where level_data: Vec>>, level_num: usize, ) -> Vec>> { - if level_num == 0 { return level_data; } @@ -208,14 +224,116 @@ where return Self::value(result, level_num - 1); } - + + fn new(data: RadarData2d>, 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>> = + 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 + FromPrimitive> MultiDimensionData for RadarData2d { -// fn levels(&self, levels: usize) -> Self { -// let a: Vec<_> = Array1::linspace(0.0, 1.0, levels) -// .into_iter() -// .map(|scale| self.resample(scale, scale, 3.0)) -// .collect(); -// } -// } +pub fn levels(data: RadarData2d>, levels: usize) -> Vec> +where + T: Num + Clone + AsPrimitive + FromPrimitive, +{ + (0..levels) + .into_iter() + .map(|level| LevelData::new(data.clone(), level)) + .collect() +} + +impl MultiDimensionData for RadarData2d +where + T: Num + Clone, + Raw: ndarray::Data + Clone + RawDataClone, +{ +} + +pub trait DataLoader { + fn load>(&self, path: P) -> Result; +} + +pub struct Npz; + +impl Npz { + #[inline] + fn load_1d( + &self, + // path: &Path, + data: &mut NpzArchive>, + name: &str, + ) -> Result, DataError> { + // let mut data = npyz::npz::NpzArchive::open(path)?; + let a = data.by_name(name)?.unwrap(); + let b: Vec = a.into_vec().unwrap(); + + Ok(Array1::from_shape_vec(b.len(), b).unwrap()) + } + #[inline] + fn load_2d( + &self, + // path: &Path, + data: &mut NpzArchive>, + name: &str, + ) -> Result, 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 = 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>) {} +} + +impl DataLoader>> for Npz +where + T: Num + Clone + Deserialize, +{ + fn load>(&self, path: P) -> Result>, DataError> { + let mut data: NpzArchive> = 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::(&mut data, "val")?, + 8 => self.load_2d::(&mut data, "val")?, + _ => {} + }, + + TypeChar::Int => {} + + _ => {} + } + } + let dim1 = self.load_1d::(&mut data, "dim1")?; + let dim2 = self.load_1d::(&mut data, "dim2")?; + let value = self.load_2d::(&mut data, "val")?; + + Ok(RadarData2d { + dim1: dim1, + dim2: dim2, + data: value, + coord_type: CoorType::LatLon, + }) + } +}