diff --git a/pyproject.toml b/pyproject.toml
index b35863e..54c9017 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "yirgacheffe"
-version = "0.8.5"
+version = "0.8.6"
 description = "Abstraction of gdal datasets for doing basic math operations"
 readme = "README.md"
 authors = [{ name = "Michael Dales", email = "mwd24@cam.ac.uk" }]
diff --git a/tests/test_operators.py b/tests/test_operators.py
index 67795cc..fe73e5c 100644
--- a/tests/test_operators.py
+++ b/tests/test_operators.py
@@ -12,13 +12,16 @@
 from yirgacheffe.operators import LayerOperation
 
 def test_add_byte_layers() -> None:
-    data1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
-    data2 = np.array([[10, 20, 30, 40], [50, 60, 70, 80]])
+    data1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]]).astype(np.uint8)
+    data2 = np.array([[10, 20, 30, 40], [50, 60, 70, 80]]).astype(np.uint8)
 
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
     layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
     result = RasterLayer.empty_raster_layer_like(layer1)
 
+    assert layer1.datatype == gdal.GDT_Byte
+    assert layer2.datatype == gdal.GDT_Byte
+
     comp = layer1 + layer2
     comp.save(result)
 
@@ -32,13 +35,16 @@ def test_add_byte_layers() -> None:
     (2, [0.0, 1.0]),
 ])
 def test_add_byte_layers_with_callback(skip, expected_steps) -> None:
-    data1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
-    data2 = np.array([[10, 20, 30, 40], [50, 60, 70, 80]])
+    data1 = np.array([[1, 2, 3, 4], [5, 6, 7, 8]]).astype(np.uint8)
+    data2 = np.array([[10, 20, 30, 40], [50, 60, 70, 80]]).astype(np.uint8)
 
     layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
     layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
     result = RasterLayer.empty_raster_layer_like(layer1)
 
+    assert layer1.datatype == gdal.GDT_Byte
+    assert layer2.datatype == gdal.GDT_Byte
+
     callback_possitions = []
 
     comp = layer1 + layer2
@@ -517,3 +523,79 @@ def test_sum_float32(monkeypatch) -> None:
             m.setattr(yirgacheffe.constants, "YSTEP", blocksize)
             actual = layer1.sum()
             assert expected == actual
+
+def test_and_byte_layers() -> None:
+    data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]]).astype(np.uint8)
+    data2 = np.array([[1, 1, 0, 0], [2, 2, 2, 2]]).astype(np.uint8)
+
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    assert layer1.datatype == gdal.GDT_Byte
+    assert layer2.datatype == gdal.GDT_Byte
+
+    comp = layer1 & layer2
+    comp.save(result)
+
+    expected = data1 & data2
+    actual = result.read_array(0, 0, 4, 2)
+
+    assert (expected == actual).all()
+
+def test_or_byte_layers() -> None:
+    data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]]).astype(np.uint8)
+    data2 = np.array([[1, 1, 0, 0], [2, 2, 2, 2]]).astype(np.uint8)
+
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    assert layer1.datatype == gdal.GDT_Byte
+    assert layer2.datatype == gdal.GDT_Byte
+
+    comp = layer1 | layer2
+    comp.save(result)
+
+    expected = data1 | data2
+    actual = result.read_array(0, 0, 4, 2)
+
+    assert (expected == actual).all()
+
+def test_and_int_layers() -> None:
+    data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]]).astype(np.int16)
+    data2 = np.array([[1, 1, 0, 0], [2, 2, 2, 2]]).astype(np.int16)
+
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    assert layer1.datatype == gdal.GDT_Int16
+    assert layer2.datatype == gdal.GDT_Int16
+
+    comp = layer1 & layer2
+    comp.save(result)
+
+    expected = data1 & data2
+    actual = result.read_array(0, 0, 4, 2)
+
+    assert (expected == actual).all()
+
+def test_or_int_layers() -> None:
+    data1 = np.array([[0, 1, 0, 2], [0, 0, 1, 1]]).astype(np.int16)
+    data2 = np.array([[1, 1, 0, 0], [2, 2, 2, 2]]).astype(np.int16)
+
+    layer1 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data1))
+    layer2 = RasterLayer(gdal_dataset_with_data((0.0, 0.0), 0.02, data2))
+    result = RasterLayer.empty_raster_layer_like(layer1)
+
+    assert layer1.datatype == gdal.GDT_Int16
+    assert layer2.datatype == gdal.GDT_Int16
+
+    comp = layer1 | layer2
+    comp.save(result)
+
+    expected = data1 | data2
+    actual = result.read_array(0, 0, 4, 2)
+
+    assert (expected == actual).all()
diff --git a/yirgacheffe/operators.py b/yirgacheffe/operators.py
index cdf7ac0..351fdef 100644
--- a/yirgacheffe/operators.py
+++ b/yirgacheffe/operators.py
@@ -59,6 +59,12 @@ def __gt__(self, other):
     def __ge__(self, other):
         return LayerOperation(self, np.ndarray.__ge__, other)
 
+    def __and__(self, other):
+        return LayerOperation(self, np.ndarray.__and__, other)
+
+    def __or__(self, other):
+        return LayerOperation(self, np.ndarray.__or__, other)
+
     def _eval(self, index, step, target_window=None):
         try:
             window = self.window