diff --git a/Cargo.toml b/Cargo.toml
index 2a78ec9..5f5aaf9 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -12,13 +12,13 @@ categories = ["science", "algorithms", "mathematics"]
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-ndarray = "0.15"
+ndarray = "0.16"
 num-traits = "0.2"
 thiserror = "1.0"
 
 [dev-dependencies]
 cargo-tarpaulin = "0.27"
-ndarray = {version = "0.15", features = ["approx-0_5", "rayon"] }
+ndarray = {version = "0.16", features = ["approx", "rayon"] }
 approx = "0.5" 
 criterion = "0.5"
 rand = "0.8"
diff --git a/benches/bench_interp1d.rs b/benches/bench_interp1d.rs
index 9ba7867..f819260 100644
--- a/benches/bench_interp1d.rs
+++ b/benches/bench_interp1d.rs
@@ -67,7 +67,7 @@ fn bench_interp1d_scalar_multithread(c: &mut Criterion) {
         })
     });
 
-    let query = query.into_shape((2500, 4)).unwrap();
+    let query = query.into_shape_with_order((2500, 4)).unwrap();
     let query_arr: Vec<_> = query.axis_iter(Axis(0)).collect();
     c.bench_function("1D scalar MT `interp_array`", |b| {
         b.iter(|| {
@@ -80,7 +80,7 @@ fn bench_interp1d_scalar_multithread(c: &mut Criterion) {
 
 fn bench_interp1d_array(c: &mut Criterion) {
     let data = Array::from_rand(500, (0.0, 1.0), 69)
-        .into_shape((100, 5))
+        .into_shape_with_order((100, 5))
         .unwrap();
     let interp = Interp1D::builder(data).build().unwrap();
     let query = Array::from_rand(10_000, (0.0, 99.0), 123);
@@ -102,7 +102,7 @@ fn bench_interp1d_array(c: &mut Criterion) {
         })
     });
 
-    let query = query.into_shape((2500, 4)).unwrap();
+    let query = query.into_shape_with_order((2500, 4)).unwrap();
     let query_arr: Vec<_> = query.axis_iter(Axis(0)).collect();
     c.bench_function("1D array `interp_array`", |b| {
         b.iter(|| {
diff --git a/benches/bench_interp1d_query_dim.rs b/benches/bench_interp1d_query_dim.rs
index 5f43396..fb00d50 100644
--- a/benches/bench_interp1d_query_dim.rs
+++ b/benches/bench_interp1d_query_dim.rs
@@ -13,7 +13,7 @@ fn bench_scalar_data_1d_query(c: &mut Criterion) {
     let interp = Interp1D::builder(data).build().unwrap();
     let query = Array::from_rand(10_000, (0.0, 99.0), 123);
 
-    let query = query.into_shape((2500, 4)).unwrap();
+    let query = query.into_shape_with_order((2500, 4)).unwrap();
     let query_arr: Vec<_> = query.axis_iter(Axis(0)).collect();
     c.bench_function("1D scalar `interp_array`", |b| {
         b.iter(|| {
@@ -38,7 +38,7 @@ fn bench_scalar_data_2d_query(c: &mut Criterion) {
     let interp = Interp1D::builder(data).build().unwrap();
     let query = Array::from_rand(10_000, (0.0, 99.0), 123);
 
-    let query = query.into_shape((625, 4, 4)).unwrap();
+    let query = query.into_shape_with_order((625, 4, 4)).unwrap();
     let query_arr: Vec<_> = query.axis_iter(Axis(0)).collect();
     c.bench_function("1D scalar `interp_array` 2D-query", |b| {
         b.iter(|| {
@@ -63,7 +63,7 @@ fn bench_scalar_data_3d_query(c: &mut Criterion) {
     let interp = Interp1D::builder(data).build().unwrap();
     let query = Array::from_rand(10_000, (0.0, 99.0), 123);
 
-    let query = query.into_shape((125, 5, 4, 4)).unwrap();
+    let query = query.into_shape_with_order((125, 5, 4, 4)).unwrap();
     let query_arr: Vec<_> = query.axis_iter(Axis(0)).collect();
     c.bench_function("1D scalar `interp_array` 3D-query", |b| {
         b.iter(|| {
diff --git a/benches/bench_interp2d.rs b/benches/bench_interp2d.rs
index fd251d9..aa300b4 100644
--- a/benches/bench_interp2d.rs
+++ b/benches/bench_interp2d.rs
@@ -11,7 +11,7 @@ mod rand_extensions;
 
 fn bench_interp2d_scalar(c: &mut Criterion) {
     let data = Array::from_rand(10_000, (0.0, 1.0), 42)
-        .into_shape((100, 100))
+        .into_shape_with_order((100, 100))
         .unwrap();
     let interp = Interp2D::builder(data).build().unwrap();
     let query_x = Array::from_rand(10_000, (0.0, 99.0), 123);
@@ -45,7 +45,7 @@ fn bench_interp2d_scalar(c: &mut Criterion) {
 
 fn bench_interp2d_scalar_multithread(c: &mut Criterion) {
     let data = Array::from_rand(10_000, (0.0, 1.0), 42)
-        .into_shape((100, 100))
+        .into_shape_with_order((100, 100))
         .unwrap();
     let interp = Interp2D::builder(data).build().unwrap();
     let query_x = Array::from_rand(10_000, (0.0, 99.0), 123);
@@ -67,8 +67,8 @@ fn bench_interp2d_scalar_multithread(c: &mut Criterion) {
         })
     });
 
-    let query_x = query_x.into_shape((2500, 4)).unwrap();
-    let query_y = query_y.into_shape((2500, 4)).unwrap();
+    let query_x = query_x.into_shape_with_order((2500, 4)).unwrap();
+    let query_y = query_y.into_shape_with_order((2500, 4)).unwrap();
     let query_arrx: Vec<_> = query_x.axis_iter(Axis(0)).collect();
     let query_arry: Vec<_> = query_y.axis_iter(Axis(0)).collect();
     c.bench_function("2D scalar MT `interp_array`", |b| {
@@ -85,7 +85,7 @@ fn bench_interp2d_scalar_multithread(c: &mut Criterion) {
 
 fn bench_interp2d_array(c: &mut Criterion) {
     let data = Array::from_rand(50_000, (0.0, 1.0), 69)
-        .into_shape((100, 100, 5))
+        .into_shape_with_order((100, 100, 5))
         .unwrap();
     let interp = Interp2D::builder(data).build().unwrap();
     let query_x = Array::from_rand(10_000, (0.0, 99.0), 123);
@@ -108,8 +108,8 @@ fn bench_interp2d_array(c: &mut Criterion) {
         })
     });
 
-    let query_x = query_x.into_shape((2500, 4)).unwrap();
-    let query_y = query_y.into_shape((2500, 4)).unwrap();
+    let query_x = query_x.into_shape_with_order((2500, 4)).unwrap();
+    let query_y = query_y.into_shape_with_order((2500, 4)).unwrap();
     let query_arrx: Vec<_> = query_x.axis_iter(Axis(0)).collect();
     let query_arry: Vec<_> = query_y.axis_iter(Axis(0)).collect();
     c.bench_function("2D array `interp_array`", |b| {
diff --git a/benches/bench_interp2d_query_dim.rs b/benches/bench_interp2d_query_dim.rs
index 95442f4..18b3682 100644
--- a/benches/bench_interp2d_query_dim.rs
+++ b/benches/bench_interp2d_query_dim.rs
@@ -8,7 +8,7 @@ mod rand_extensions;
 
 fn setup() -> (Interp2DScalar<f64, Bilinear>, Array1<f64>, Array1<f64>) {
     let data = Array::from_rand(10_000, (0.0, 1.0), 42)
-        .into_shape((100, 100))
+        .into_shape_with_order((100, 100))
         .unwrap();
     let interp = Interp2D::builder(data).build().unwrap();
     let query_x = Array::from_rand(10_000, (0.0, 99.0), 123);
@@ -18,8 +18,8 @@ fn setup() -> (Interp2DScalar<f64, Bilinear>, Array1<f64>, Array1<f64>) {
 
 fn bench_scalar_data_1d_query(c: &mut Criterion) {
     let (interp, query_x, query_y) = black_box(setup());
-    let query_x = query_x.into_shape((2500, 4)).unwrap();
-    let query_y = query_y.into_shape((2500, 4)).unwrap();
+    let query_x = query_x.into_shape_with_order((2500, 4)).unwrap();
+    let query_y = query_y.into_shape_with_order((2500, 4)).unwrap();
     let query_arrx: Vec<_> = query_x.axis_iter(Axis(0)).collect();
     let query_arry: Vec<_> = query_y.axis_iter(Axis(0)).collect();
 
@@ -43,8 +43,8 @@ fn bench_scalar_data_1d_query(c: &mut Criterion) {
 
 fn bench_scalar_data_2d_query(c: &mut Criterion) {
     let (interp, query_x, query_y) = black_box(setup());
-    let query_x = query_x.into_shape((625, 4, 4)).unwrap();
-    let query_y = query_y.into_shape((625, 4, 4)).unwrap();
+    let query_x = query_x.into_shape_with_order((625, 4, 4)).unwrap();
+    let query_y = query_y.into_shape_with_order((625, 4, 4)).unwrap();
     let query_arrx: Vec<_> = query_x.axis_iter(Axis(0)).collect();
     let query_arry: Vec<_> = query_y.axis_iter(Axis(0)).collect();
 
@@ -68,8 +68,8 @@ fn bench_scalar_data_2d_query(c: &mut Criterion) {
 
 fn bench_scalar_data_3d_query(c: &mut Criterion) {
     let (interp, query_x, query_y) = black_box(setup());
-    let query_x = query_x.into_shape((125, 5, 4, 4)).unwrap();
-    let query_y = query_y.into_shape((125, 5, 4, 4)).unwrap();
+    let query_x = query_x.into_shape_with_order((125, 5, 4, 4)).unwrap();
+    let query_y = query_y.into_shape_with_order((125, 5, 4, 4)).unwrap();
     let query_arrx: Vec<_> = query_x.axis_iter(Axis(0)).collect();
     let query_arry: Vec<_> = query_y.axis_iter(Axis(0)).collect();
 
diff --git a/src/interp1d/mod.rs b/src/interp1d/mod.rs
index ac45d7a..3931b03 100644
--- a/src/interp1d/mod.rs
+++ b/src/interp1d/mod.rs
@@ -308,14 +308,15 @@ where
                     }
                 });
 
-            let subview = match subview.into_shape(self.data.raw_dim().remove_axis(Axis(0))) {
-                Ok(view) => view,
-                Err(err) => {
-                    let expect = self.get_buffer_shape(xs.raw_dim()).into_pattern();
-                    let got = buffer.dim();
-                    panic!("{err} expected: {expect:?}, got: {got:?}")
-                }
-            };
+            let subview =
+                match subview.into_shape_with_order(self.data.raw_dim().remove_axis(Axis(0))) {
+                    Ok(view) => view,
+                    Err(err) => {
+                        let expect = self.get_buffer_shape(xs.raw_dim()).into_pattern();
+                        let got = buffer.dim();
+                        panic!("{err} expected: {expect:?}, got: {got:?}")
+                    }
+                };
 
             self.strategy.interp_into(self, subview, x)?;
         }
@@ -359,7 +360,11 @@ where
     ///  - `x` is stricktly monotonic rising
     ///  - `data.shape()[0] == x.len()`
     ///  - the `strategy` is porperly initialized with the data
-    pub unsafe fn new_unchecked(x: ArrayBase<Sx, Ix1>, data: ArrayBase<Sd, D>, strategy: Strat) -> Self {
+    pub unsafe fn new_unchecked(
+        x: ArrayBase<Sx, Ix1>,
+        data: ArrayBase<Sd, D>,
+        strategy: Strat,
+    ) -> Self {
         Interp1D { x, data, strategy }
     }
 
@@ -498,7 +503,7 @@ mod tests {
     macro_rules! get_interp {
         ($dim:expr, $shape:expr) => {{
             let arr = rand_arr(4usize.pow($dim), (0.0, 1.0), 64)
-                .into_shape($shape)
+                .into_shape_with_order($shape)
                 .unwrap();
             Interp1D::builder(arr).build().unwrap()
         }};
@@ -565,7 +570,7 @@ mod tests {
     #[should_panic(expected = "expected: [2], got: [1]")] // this is not really a good message
     fn interp1d_2d_array_into_too_small1() {
         let arr = rand_arr((4usize).pow(2), (0.0, 1.0), 64)
-            .into_shape((4, 4))
+            .into_shape_with_order((4, 4))
             .unwrap();
         let interp = Interp1D::builder(arr).build().unwrap();
         let mut buf = Array::zeros((1, 4));
@@ -576,7 +581,7 @@ mod tests {
     #[should_panic]
     fn interp1d_2d_array_into_too_small2() {
         let arr = rand_arr((4usize).pow(2), (0.0, 1.0), 64)
-            .into_shape((4, 4))
+            .into_shape_with_order((4, 4))
             .unwrap();
         let interp = Interp1D::builder(arr).build().unwrap();
         let mut buf = Array::zeros((2, 3));
@@ -587,7 +592,7 @@ mod tests {
     #[should_panic]
     fn interp1d_2d_array_into_too_big1() {
         let arr = rand_arr((4usize).pow(2), (0.0, 1.0), 64)
-            .into_shape((4, 4))
+            .into_shape_with_order((4, 4))
             .unwrap();
         let interp = Interp1D::builder(arr).build().unwrap();
         let mut buf = Array::zeros((3, 4));
@@ -598,7 +603,7 @@ mod tests {
     #[should_panic]
     fn interp1d_2d_array_into_too_big2() {
         let arr = rand_arr((4usize).pow(2), (0.0, 1.0), 64)
-            .into_shape((4, 4))
+            .into_shape_with_order((4, 4))
             .unwrap();
         let interp = Interp1D::builder(arr).build().unwrap();
         let mut buf = Array::zeros((2, 5));
diff --git a/src/interp2d/mod.rs b/src/interp2d/mod.rs
index e516caf..87c0b87 100644
--- a/src/interp2d/mod.rs
+++ b/src/interp2d/mod.rs
@@ -265,7 +265,7 @@ where
                     }
                 });
 
-            let subview = match subview.into_shape(
+            let subview = match subview.into_shape_with_order(
                 self.data
                     .raw_dim()
                     .remove_axis(Axis(0))
@@ -543,7 +543,7 @@ mod tests {
             #[test]
             fn $name() {
                 let arr = rand_arr(4usize.pow($dim), (0.0, 1.0), 64)
-                    .into_shape($shape)
+                    .into_shape_with_order($shape)
                     .unwrap();
                 let interp = Interp2D::builder(arr).build().unwrap();
                 let res = interp.interp(2.2, 2.2).unwrap();
@@ -578,7 +578,7 @@ mod tests {
     #[test]
     fn interp2d_2d_scalar() {
         let arr = rand_arr(4usize.pow(2), (0.0, 1.0), 64)
-            .into_shape((4, 4))
+            .into_shape_with_order((4, 4))
             .unwrap();
         let _res: f64 = Interp2D::builder(arr) // typecheck f64 as return type
             .build()
diff --git a/tests/interp2d.rs b/tests/interp2d.rs
index 4da56e8..6efd5bc 100644
--- a/tests/interp2d.rs
+++ b/tests/interp2d.rs
@@ -83,17 +83,19 @@ fn extrapolate() {
 
 #[test]
 fn interpolate_array() {
-    let data = Array::linspace(0.0, 8.0, 9).into_shape((3, 3)).unwrap();
+    let data = Array::linspace(0.0, 8.0, 9)
+        .into_shape_with_order((3, 3))
+        .unwrap();
     let x = array![1.0, 2.0, 3.0];
     let y = array![4.0, 5.0, 6.0];
     let resolution = 11usize;
     let qx = Array::linspace(1.0, 3.0, resolution);
     let qy = Array::linspace(4.0, 6.0, resolution);
     let qx = Array::from_iter(qx.into_iter().flat_map(|x| repeat(x).take(resolution)))
-        .into_shape((resolution, resolution))
+        .into_shape_with_order((resolution, resolution))
         .unwrap();
     let qy = Array::from_iter(repeat(qy).take(resolution).flatten())
-        .into_shape((resolution, resolution))
+        .into_shape_with_order((resolution, resolution))
         .unwrap();
 
     let interp = Interp2D::builder(data).x(x).y(y).build().unwrap();
@@ -247,7 +249,7 @@ fn interp_nd_data() {
         .flatten()
         .flatten(),
     )
-    .into_shape((2, 2, 2, 2))
+    .into_shape_with_order((2, 2, 2, 2))
     .unwrap();
 
     let interp = Interp2DBuilder::new(data).build().unwrap();
@@ -265,7 +267,9 @@ fn interp_nd_data() {
 #[test]
 #[should_panic(expected = "`xs.shape()` and `ys.shape()` do not match")]
 fn interp_array_with_unmatched_axis() {
-    let data = Array::linspace(0.0, 8.0, 9).into_shape((3, 3)).unwrap();
+    let data = Array::linspace(0.0, 8.0, 9)
+        .into_shape_with_order((3, 3))
+        .unwrap();
     let qx = array![0.0, 1.0];
     let qy = array![0.0, 1.0, 2.0];
     let interp = Interp2D::builder(data).build().unwrap();