f.data[0:0] != np.array(f.data)[0:0]
Reproducer

```
In [8]: from devito import Grid, Function                                                                                                                                                                          

In [9]: grid = Grid(shape=(4, 4))                                                                                                                                                                                  

In [10]: f = Function(name='f', grid=grid)                                                                                                                                                                         

In [11]: f.data                                                                                                                                                                                                    
Allocating memory for f(6, 6)
Out[11]: 
Data([[0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.]], dtype=float32)

In [12]: f.data[0:0]                                                                                                                                                                                               
Out[12]: 
Data([[0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.]], dtype=float32)

In [13]: f.data[slice(0, 0)]                                                                                                                                                                                       
Out[13]: 
Data([[0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.]], dtype=float32)

In [14]: np.array(f.data)                                                                                                                                                                                          
Out[14]: 
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]], dtype=float32)

In [15]: np.array(f.data)[0:0]                                                                                                                                                                                     
Out[15]: array([], shape=(0, 4), dtype=float32)

In [16]: np.array(f.data)[slice(0, 0)]                                                                                                                                                                             
Out[16]: array([], shape=(0, 4), dtype=float32)
```

issue cropping up at In [12], as f.data[0:0] should return an empty array
also note that the problem only seems to be 0:0 :

```
In [19]: f.data[0:0]                                                                                                                                                                                               
Out[19]: 
Data([[0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.]], dtype=float32)

In [20]: f.data[1:1]                                                                                                                                                                                               
Out[20]: Data([], shape=(0, 4), dtype=float32)

In [21]: f.data[3:3]                                                                                                                                                                                               
Out[21]: Data([], shape=(0, 4), dtype=float32)

In [22]: f.data[-1:-1]                                                                                                                                                                                             
Out[22]: Data([], shape=(0, 4), dtype=float32)

In [23]: f.data[:, 0:0]                                                                                                                                                                                            
Out[23]: 
Data([[0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.],
      [0., 0., 0., 0.]], dtype=float32)

In [24]: f.data[:, -1:-1]                                                                                                                                                                                          
Out[24]: Data([], shape=(4, 0), dtype=float32)

In [25]: f.data[:, 2:2]                                                                                                                                                                                            
Out[25]: Data([], shape=(4, 0), dtype=float32)
15:41
(in any dimension)
```
diff --git a/devito/builtins/initializers.py b/devito/builtins/initializers.py
index 1738e0e7a..b829863fe 100644
--- a/devito/builtins/initializers.py
+++ b/devito/builtins/initializers.py
@@ -283,9 +283,7 @@ def initialize_function(function, data, nbl, mapper=None, mode='constant',
             pad_outhalo(function)
         return
 
-    nbl = nbl_to_padsize(nbl, function.ndim)
-    slices = tuple([slice(nl, -nr) for _, (nl, nr) in zip(range(function.grid.dim),
-                                                          as_tuple(nbl))])
+    nbl, slices = nbl_to_padsize(nbl, function.ndim)
     if isinstance(data, dv.Function):
         function.data[slices] = data.data[:]
     else:
diff --git a/devito/builtins/utils.py b/devito/builtins/utils.py
index 09a46444e..64f54d24b 100644
--- a/devito/builtins/utils.py
+++ b/devito/builtins/utils.py
@@ -50,18 +50,27 @@ def nbl_to_padsize(nbl, ndim):
     nb_total = as_tuple(nbl)
     if len(nb_total) == 1 and len(nb_total) < ndim:
         nb_pad = ndim*((nb_total[0], nb_total[0]),)
+        slices = ndim*(slice(nb_total[0], None, 1),) \\
+            if nb_total[0] == 0 else ndim*(slice(nb_total[0], -nb_total[0], 1),)
     elif len(nb_total) == ndim:
         nb_pad = []
+        slices = []
         for nb_dim in nb_total:
             if len(as_tuple(nb_dim)) == 1:
                 nb_pad.append((nb_dim, nb_dim))
+                s = slice(nb_dim, None, 1) if nb_dim == 0 \\
+                    else slice(nb_dim, -nb_dim, 1)
+                slices.append(s)
             else:
                 assert len(nb_dim) == 2
                 nb_pad.append(nb_dim)
+                s = slice(nb_dim[0], None, 1) if nb_dim[1] == 0 \\
+                    else slice(nb_dim[0], -nb_dim[1], 1)
+                slices.append(s)
     else:
         raise ValueError("`nbl` must be an integer or tuple of (tuple of) integers"
                          "of length `function.ndim`.")
-    return tuple(nb_pad)
+    return tuple(nb_pad), tuple(slices)
 
 
 def pad_outhalo(function):
diff --git a/devito/data/data.py b/devito/data/data.py
index b91e23d07..ad4506872 100644
--- a/devito/data/data.py
+++ b/devito/data/data.py
@@ -415,8 +415,7 @@ class Data(np.ndarray):
                 # Need to wrap index based on modulo
                 v = index_apply_modulo(i, s)
             elif self._is_distributed is True and dec is not None:
-                # Need to convert the user-provided global indices into local indices.
-                # Obviously this will have no effect if MPI is not used
+                # Convert the user-provided global indices into local indices.
                 try:
                     v = convert_index(i, dec, mode='glb_to_loc')
                 except TypeError:
diff --git a/devito/data/decomposition.py b/devito/data/decomposition.py
index c157feed9..2bb3faf79 100644
--- a/devito/data/decomposition.py
+++ b/devito/data/decomposition.py
@@ -239,7 +239,12 @@ class Decomposition(tuple):
                 elif isinstance(glb_idx, slice):
                     if self.loc_empty:
                         return slice(-1, -3)
-                    if glb_idx.step >= 0:
+                    if glb_idx.step >= 0 and glb_idx.stop == self.glb_min:
+                        glb_idx_min = self.glb_min if glb_idx.start is None \\
+                            else glb_idx.start
+                        glb_idx_max = self.glb_min
+                        retfunc = lambda a, b: slice(a, b, glb_idx.step)
+                    elif glb_idx.step >= 0:
                         glb_idx_min = self.glb_min if glb_idx.start is None \\
                             else glb_idx.start
                         glb_idx_max = self.glb_max if glb_idx.stop is None \\
diff --git a/tests/test_data.py b/tests/test_data.py
index d3222da0d..7bc96e35a 100644
--- a/tests/test_data.py
+++ b/tests/test_data.py
@@ -1034,6 +1034,21 @@ class TestDataDistributed(object):
         expected2 = np.full(result2.shape, 1)
         assert(np.all(result2 == expected2))
 
+    def test_empty_slicing(self):
+        grid = Grid(shape=(2, 2), extent=(1, 1))
+        f = Function(name='f', grid=grid)
+        g = TimeFunction(name='g', grid=grid)
+        assert(f.data[0:0].shape == (0, 2))
+        assert(f.data[0:0, 0:0].shape == (0, 0))
+        assert(g.data[0:0].shape == (0, 2, 2))
+        assert(g.data[0:0, 0:0].shape == (0, 0, 2))
+        assert(g.data[0:0, 0:0, 0:0].shape == (0, 0, 0))
+        assert(f.data[1:1].shape == (0, 2))
+        assert(f.data[0:0, 1:1].shape == (0, 0))
+        assert(g.data[1:1].shape == (0, 2, 2))
+        assert(g.data[1:1, 0:0].shape == (0, 0, 2))
+        assert(g.data[1:1, 0:0, 1:1].shape == (0, 0, 0))
+
     @pytest.mark.parallel(mode=4)
     def test_neg_start_stop(self):
         grid0 = Grid(shape=(8, 8))
['devito/data/data.py', 'devito/data/decomposition.py', 'tests/test_data.py', 'devito/builtins/utils.py', 'devito/builtins/initializers.py']