diff --git a/karta/raster/read.py b/karta/raster/read.py
index 3f0d9a7..8f75018 100644
--- a/karta/raster/read.py
+++ b/karta/raster/read.py
@@ -75,13 +75,19 @@ class of band used by returned grid (default karta.band.CompressedBand)
         will have bands of type karta.raster._gdal.GdalFileBand
     """
     in_memory = kw.get("in_memory", True)
+    if len(fnms) == 0:
+        raise ValueError("at least one filename must be provided")
     bands = []
     hdrs = []
-    for fnm in fnms:
+    for i, fnm in enumerate(fnms):
         _b, _h = _gdal.read(fnm, in_memory, 1, **kw)
-        if len(hdrs) != 0:
-            if _h != hdrs[-1]:
-                raise GridError("geotransform in {0} not equivalent to a previous GeoTiff".format(fnm))
+        if (i != 0) and (len(hdrs) != 0):
+            for k, v in _h.items():
+                if (isinstance(v, float)
+                        and np.isnan(v) and not np.isnan(hdrs[-1].get(k, 0))):
+                    if v != hdrs[-1].get(k, np.nan):
+                        raise GridError("geotransform in {0} not equivalent to "
+                                        "a previous GeoTiff".format(fnm))
         bands.append(_b[0])
         hdrs.append(_h)
 
diff --git a/tests/geotiff_tests.py b/tests/geotiff_tests.py
index 8f2d7bd..130a958 100644
--- a/tests/geotiff_tests.py
+++ b/tests/geotiff_tests.py
@@ -21,7 +21,7 @@ def test_numpy_type_coercion(self):
         self.assertEqual(_gdal.numpy_dtype(11), np.complex64)
         return
 
-    def test_io(self):
+    def test_write_read(self):
         # try writing a file, then read it back in and verify that it matches
         v = peaks(500)[:100,:]
         utm7 = karta.crs.ProjectedCRS("+proj=utm +zone=7 +north +datum=WGS84",
@@ -39,7 +39,7 @@ def test_io(self):
         self.assertTrue(np.all(g[:,:] == gnew[:,:]))
         return
 
-    def test_io_virtual(self):
+    def test_write_read_disk(self):
         # try writing a file, then open it without loading into memory and verify
         v = peaks(500)[:100,:]
         utm7 = karta.crs.ProjectedCRS("+proj=utm +zone=7 +north +datum=WGS84",
@@ -65,10 +65,39 @@ def test_write_compress(self):
         g = karta.RegularGrid([15.0, 15.0, 30.0, 30.0, 0.0, 0.0], v, crs=utm7)
 
         fpath = os.path.join(TMPDATA, "test.tif")
-        g.to_gtiff(fpath, compress="LZW")
-        g.to_gtiff(fpath, compress="PACKBITS")
+        g.to_geotiff(fpath, compress="LZW")
+        self.assertTrue(os.path.isfile(fpath))
+        os.remove(fpath)
+        g.to_geotiff(fpath, compress="PACKBITS")
+        self.assertTrue(os.path.isfile(fpath))
         return
 
+    def test_read_as_bands(self):
+        # write several files and then read as a single multiband grid
+        v = peaks(500)[:100,:]
+        utm7 = karta.crs.ProjectedCRS("+proj=utm +zone=7 +north +datum=WGS84",
+                                      "UTM 7N (WGS 84)")
+        g1 = karta.RegularGrid([15.0, 15.0, 30.0, 30.0, 0.0, 0.0], v, crs=utm7)
+        g2 = karta.RegularGrid([15.0, 15.0, 30.0, 30.0, 0.0, 0.0], v**2, crs=utm7)
+        g3 = karta.RegularGrid([15.0, 15.0, 30.0, 30.0, 0.0, 0.0], v+2, crs=utm7)
+        g4 = karta.RegularGrid([15.0, 15.0, 30.0, 30.0, 0.0, 0.0], v*2, crs=utm7)
+
+        paths = []
+        for i, g in enumerate((g1, g2, g3, g4)):
+            fpath = os.path.join(TMPDATA, "test{0}.tif".format(i))
+            g.to_geotiff(fpath, compress=None)
+            paths.append(fpath)
+
+        gnew = karta.from_geotiffs(*paths)
+        self.assertTrue("+proj=utm" in gnew.crs.get_proj4())
+        self.assertTrue("+zone=7" in gnew.crs.get_proj4())
+        self.assertEqual(g.transform, gnew.transform)
+        self.assertEqual(g.values.dtype, gnew.values.dtype)
+        self.assertEqual(gnew.size, (100, 500))
+        self.assertEqual(gnew.nbands, 4)
+        return
+
+
 class GdalVirtualArrayTests(unittest.TestCase):
 
     def setUp(self):
@@ -82,6 +111,7 @@ def setUp(self):
         self.grid = karta.read_gtiff(fpath, in_memory=False)
 
     def test_slicing_virtual(self):
+        """ make sure that slicing a disk-based array works """
         self.grid[5:10, 7:15]
         self.grid[5:10:2, 7:15]
         self.grid[10:5:-1, 7:15]
@@ -96,9 +126,11 @@ def test_slicing_virtual(self):
         return
 
     def test_iteration_virtual(self):
+        i = 0
         for row in self.grid.values:
-            pass
-
+            i += 1
+        self.assertEqual(i, 100)
+        return
 
 def peaks(n=49):
     """ 2d peaks function of MATLAB logo fame. """