Repository: DHI/mikeio

Issue:
isel() always returns a copy of the data
If you want to change a subset of the values in a DataArray, you may try something like: 

     da.isel(time=0).values[:] = 7.0 

or 

    da[0].values[:] = 7.0 

This will unfortunately not work as isel() (and getitem) returns a *copy* of the data even when it does not need to. The general rule in NumPy is that simple indexing (single index or slice) will return a view, fancy indexing (e.g. a list of ids or boolean index) will return a copy. MIKE IO should follow the same pattern. 

Until this is fixed, you should be able to reformulate your code like this: 

    da.values[0,:] = 7.0 
 
Thanks to @msaes88 for reporting! 

Discussion:

Patch:
diff --git a/docs/dataarray.md b/docs/dataarray.md
index ad4b0073..41a223ec 100644
--- a/docs/dataarray.md
+++ b/docs/dataarray.md
@@ -60,6 +60,33 @@ geometry: GeometryPoint2D(x=607002.7094112666, y=6906734.833048992)
 values: [0.4591, 0.8078, ..., -0.6311]
 ```
 
+## Modifying values
+
+You can modify the values of a DataArray by changing its `values`: 
+
+```python
+>>> da.values[0, 3] = 5.0
+```
+
+If you wish to change values of a subset of the DataArray you should be aware of the difference between a _view_ and a _copy_ of the data. Similar to NumPy, MIKE IO selection method will return a _view_ of the data when using single index and slices, but a _copy_ of the data using fancy indexing (a list of indicies or boolean indexing). Note that prior to release 1.3, MIKE IO would always return a copy. 
+
+It is recommended to change the values using `values` property directly on the original DataArray (like above), but it is also possible to change the values of the original DataArray by working on a subset DataArray if it is selected with single index or slice as explained above. 
+
+```python
+>>> da_sub = da.isel(time=0)
+>>> da_sub.values[:] = 5.0    # will change da
+```
+
+Fancy indexing will return a _copy_ and therefore not change the original:
+
+```python
+>>> da_sub = da.isel(time=[0,1,2])
+>>> da_sub.values[:] = 5.0    # will NOT change da
+```
+
+
+
+
 ## Plotting
 
 The plotting of a DataArray is context-aware meaning that plotting behaviour depends on the geometry of the DataArray being plotted. 
diff --git a/mikeio/data_utils.py b/mikeio/data_utils.py
index e1e48241..d615f3c1 100644
--- a/mikeio/data_utils.py
+++ b/mikeio/data_utils.py
@@ -53,12 +53,19 @@ class DataUtilsMixin:
                 )
                 steps = list(range(s.start, s.stop))
             except TypeError:
-                steps = list(range(*steps.indices(len(time))))
+                pass
+                #steps = list(range(*steps.indices(len(time))))
         elif isinstance(steps, int):
             steps = [steps]
 
         return steps
 
+    @staticmethod
+    def _n_selected_timesteps(time, k):
+        if isinstance(k, slice):
+            k = list(range(*k.indices(len(time))))
+        return len(k)
+            
     @staticmethod
     def _is_boolean_mask(x) -> bool:
         if hasattr(x, "dtype"):  # isinstance(x, (np.ndarray, DataArray)):
diff --git a/mikeio/dataarray.py b/mikeio/dataarray.py
index 3e30d0e9..818cf721 100644
--- a/mikeio/dataarray.py
+++ b/mikeio/dataarray.py
@@ -1013,7 +1013,10 @@ class DataArray(DataUtilsMixin, TimeSeries):
 
     @values.setter
     def values(self, value):
-        if value.shape != self._values.shape:
+        if np.isscalar(self._values):
+            if not np.isscalar(value):
+                raise ValueError("Shape of new data is wrong (should be scalar)")
+        elif value.shape != self._values.shape:
             raise ValueError("Shape of new data is wrong")
 
         self._values = value
@@ -1095,7 +1098,7 @@ class DataArray(DataUtilsMixin, TimeSeries):
                 if dims[j] == "time":
                     # getitem accepts fancy indexing only for time
                     k = self._get_time_idx_list(self.time, k)
-                    if len(k) == 0:
+                    if self._n_selected_timesteps(self.time, k) == 0:
                         raise IndexError("No timesteps found!")
                 da = da.isel(k, axis=dims[j])
         return da
@@ -1137,6 +1140,10 @@ class DataArray(DataUtilsMixin, TimeSeries):
         """Return a new DataArray whose data is given by
         integer indexing along the specified dimension(s).
 
+        Note that the data will be a _view_ of the original data 
+        if possible (single index or slice), otherwise a copy (fancy indexing) 
+        following NumPy convention. 
+
         The spatial parameters available depend on the dims
         (i.e. geometry) of the DataArray:
 
@@ -1238,7 +1245,9 @@ class DataArray(DataUtilsMixin, TimeSeries):
 
         axis = self._parse_axis(self.shape, self.dims, axis)
 
+        idx_slice = None
         if isinstance(idx, slice):
+            idx_slice = idx
             idx = list(range(*idx.indices(self.shape[axis])))
         if idx is None or (not np.isscalar(idx) and len(idx) == 0):
             return None
@@ -1264,12 +1273,22 @@ class DataArray(DataUtilsMixin, TimeSeries):
                 )
                 zn = self._zn[:, node_ids]
 
+        # reduce dims only if singleton idx
+        dims = tuple([d for i, d in enumerate(self.dims) if i != axis]) if single_index else self.dims
         if single_index:
-            # reduce dims only if singleton idx
-            dims = tuple([d for i, d in enumerate(self.dims) if i != axis])
-            dat = np.take(self.values, int(idx), axis=axis)
+            idx = int(idx)
+        elif idx_slice is not None:
+            idx = idx_slice
+
+        if axis == 0:
+            dat = self.values[idx]
+        elif axis == 1:
+            dat = self.values[:,idx]
+        elif axis == 2:
+            dat = self.values[:,:,idx]
+        elif axis == 3:
+            dat = self.values[:,:,:,idx]
         else:
-            dims = self.dims
             dat = np.take(self.values, idx, axis=axis)
 
         return DataArray(
diff --git a/mikeio/dataset.py b/mikeio/dataset.py
index 39558816..7e2ba67b 100644
--- a/mikeio/dataset.py
+++ b/mikeio/dataset.py
@@ -698,7 +698,7 @@ class Dataset(DataUtilsMixin, TimeSeries, collections.abc.MutableMapping):
             key = pd.DatetimeIndex(key)
         if isinstance(key, pd.DatetimeIndex) or self._is_key_time(key):
             time_steps = self._get_time_idx_list(self.time, key)
-            if len(time_steps) == 0:
+            if self._n_selected_timesteps(self.time, time_steps) == 0:
                 raise IndexError("No timesteps found!")
             return self.isel(time_steps, axis=0)
         if isinstance(key, slice):


Test patch:
diff --git a/tests/test_dataarray.py b/tests/test_dataarray.py
index 7dae466e..56b33dde 100644
--- a/tests/test_dataarray.py
+++ b/tests/test_dataarray.py
@@ -839,6 +839,98 @@ def test_modify_values(da1):
     da1.values = np.zeros_like(da1.values) + 2.0
 
 
+def test_modify_values_1d(da1):
+    assert da1.values[4] == 14.0
+
+    # selecting a slice will return a view. The original is changed.
+    da1.isel(slice(4,6)).values[0] = 13.0 
+    assert da1.values[4] == 13.0
+
+    # __getitem__ uses isel()
+    da1[4:6].values[0] = 12.0 
+    assert da1.values[4] == 12.0
+
+    # values is scalar, therefore copy by definition. Original is not changed.
+    da1.isel(4).values = 11.0    
+    assert da1.values[4] != 11.0  
+
+    # fancy indexing will return copy! Original is *not* changed.
+    da1.isel([0,4,7]).values[1] = 10.0  
+    assert da1.values[4] != 10.0  
+
+
+def test_modify_values_2d_all(da2):
+    assert da2.shape == (10,7)
+    assert da2.values[2,5] == 0.1
+
+    da2 += 0.1
+    assert da2.values[2,5] == 0.2
+
+    vals = 0.3*np.ones(da2.shape)
+    da2.values = vals
+    assert da2.values[2,5] == 0.3
+
+
+def test_modify_values_2d_idx(da2):
+    assert da2.shape == (10,7)
+    assert da2.values[2,5] == 0.1
+
+    # selecting a single index will return a view. The original is changed.
+    da2.isel(time=2).values[5] = 0.2
+    assert da2.values[2,5] == 0.2
+
+    da2.isel(x=5).values[2] = 0.3
+    assert da2.values[2,5] == 0.3
+
+    da2.values[2,5] = 0.4
+    assert da2.values[2,5] == 0.4
+
+    # __getitem__ uses isel()
+    da2[2].values[5] = 0.5
+    assert da2.values[2,5] == 0.5
+
+    da2[:,5].values[2] = 0.6
+    assert da2.values[2,5] == 0.6
+
+
+def test_modify_values_2d_slice(da2):
+    assert da2.shape == (10,7)
+    assert da2.values[2,5] == 0.1
+
+    # selecting a slice will return a view. The original is changed.
+    da2.isel(time=slice(2,6)).values[0,5] = 0.4
+    assert da2.values[2,5] == 0.4
+
+    da2.isel(x=slice(5,7)).values[2,0] = 0.5
+    assert da2.values[2,5] == 0.5
+
+    # __getitem__ uses isel()
+    da2[2:5].values[0,5] = 0.6
+    assert da2.values[2,5] == 0.6
+
+    da2[:,5:7].values[2,0] = 0.7
+    assert da2.values[2,5] == 0.7
+
+
+def test_modify_values_2d_fancy(da2):
+    assert da2.shape == (10,7)
+    assert da2.values[2,5] == 0.1
+
+    # fancy indexing will return a *copy*. The original is NOT changed.
+    da2.isel(time=[2,3,4,5]).values[0,5] = 0.4
+    assert da2.values[2,5] != 0.4
+
+    da2.isel(x=[5,6]).values[2,0] = 0.5
+    assert da2.values[2,5] != 0.5
+
+    # __getitem__ uses isel()
+    da2[[2,3,4,5]].values[0,5] = 0.6
+    assert da2.values[2,5] != 0.6
+
+    da2[:,[5,6]].values[2,0] = 0.7
+    assert da2.values[2,5] != 0.7
+
+
 def test_add_scalar(da1):
     da2 = da1 + 10.0
     assert isinstance(da2, mikeio.DataArray)
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index c0edc778..fd8e8239 100644
--- a/tests/test_dataset.py
+++ b/tests/test_dataset.py
@@ -273,7 +273,7 @@ def test_select_subset_isel_axis_out_of_range_error(ds2):
     # After subsetting there is only one dimension
     assert len(dss.shape) == 1
 
-    with pytest.raises(ValueError):
+    with pytest.raises(IndexError):
         dss.isel(idx=0, axis=1)