diff --git a/paraview-dns-data-reader.py b/paraview-dns-data-reader.py
index 010f7ca..edb809e 100644
--- a/paraview-dns-data-reader.py
+++ b/paraview-dns-data-reader.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+import vtk
 from vtkmodules.vtkCommonDataModel import vtkDataSet
 from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase
 from vtkmodules.numpy_interface import dataset_adapter as dsa
@@ -17,7 +18,7 @@ class PythonDnsDataReader(VTKPythonAlgorithmBase):
         VTKPythonAlgorithmBase.__init__(self,
                 nInputPorts=0,
                 nOutputPorts=1,
-                outputType='vtkRectilinearGrid')
+                outputType='vtkImageData')
         self._filename = None
 
     @smproperty.stringvector(name="FileName")
@@ -31,11 +32,26 @@ class PythonDnsDataReader(VTKPythonAlgorithmBase):
             self._timesteps = None
             self.Modified()
 
+    def RequestInformation(self, request, inInfoVec, outInfoVec):
+        executive = self.GetExecutive()
+        outInfo = outInfoVec.GetInformationObject(0)
+
+        print("RequestInformation: ", self._filename)
+
+        filename = self._filename
+        t, nx, ny, nz, U, V, W, Y0, Y1 = self._read_data(filename, True)
+
+        dx = 4*np.pi/nx
+        dy = 2*np.pi/ny
+        dz = 2*np.pi/nz
+
+        outInfo.Set(executive.WHOLE_EXTENT(), 0, nx-1, 0, ny-1, 0, nz-1)
+        outInfo.Set(vtk.vtkDataObject.SPACING(), dx, dy, dz)
+
+        return 1
 
     def RequestData(self, request, inInfo, outInfo):
-        import vtk
-        from vtkmodules.vtkCommonDataModel import vtkRectilinearGrid
-        output = vtkRectilinearGrid.GetData(outInfo, 0)
+        output = dsa.WrapDataObject(self.GetOutputData(outInfo, 0))
 
         if self._filename is None:
             # Note, exceptions are totally fine!
@@ -51,60 +67,12 @@ class PythonDnsDataReader(VTKPythonAlgorithmBase):
         w = W.ravel()
         y1 = Y1.ravel()
 
-        x = np.arange(nx) * 4*np.pi/nx
-        y = np.arange(ny) * 2*np.pi/ny
-        z = np.arange(nz) * 2*np.pi/nz
+        dims = [nx, ny, nz]
+        assert y1.shape[0] == dims[0]*dims[1]*dims[2], "dimension mismatch"
 
-        # Create a rectilinear grid by defining three arrays specifying the
-        # coordinates in the x-y-z directions.
-        xCoords = vtk.vtkFloatArray()
-        xCoords.SetNumberOfTuples(len(x))
-        for i, xi in enumerate(x):
-            xCoords.SetTuple1(i, xi)
-
-        yCoords = vtk.vtkFloatArray()
-        yCoords.SetNumberOfTuples(len(y))
-        for i, yi in enumerate(y):
-            yCoords.SetTuple1(i, yi)
-
-        zCoords = vtk.vtkFloatArray()
-        zCoords.SetNumberOfTuples(len(z))
-        for i, zi in enumerate(z):
-            zCoords.SetTuple1(i, zi)
-
-        # The coordinates are assigned to the rectilinear grid. Make sure that
-        # the number of values in each of the XCoordinates, YCoordinates,
-        # and ZCoordinates is equal to what is defined in SetDimensions().
-        #
-        output.SetDimensions(len(x), len(y), len(z))
-        output.SetXCoordinates(xCoords)
-        output.SetYCoordinates(yCoords)
-        output.SetZCoordinates(zCoords)
-
-        numPoints = output.GetNumberOfPoints()
-
-        velocity = vtk.vtkFloatArray();
-        velocity.SetNumberOfComponents(3);
-        velocity.SetNumberOfTuples(numPoints);
-        velocity.SetName("U");
-
-        for i in range(0, numPoints):
-            velocity.SetTuple3(i, u[i], v[i], w[i]);
-
-        output.GetPointData().AddArray(velocity)
-        del velocity
-
-        c = vtk.vtkFloatArray();
-        c.SetNumberOfTuples(numPoints);
-        c.SetName("c");
-
-        for i in range(0, numPoints):
-            c.SetTuple1(i, 1.0-y1[i]);
-
-        output.GetPointData().AddArray(c)
-        del c
-
-        del t, nx, ny, nz, u, v, w, Y0, Y1, x, y, z
+        output.SetExtent(0, dims[0]-1, 0, dims[1]-1, 0, dims[2]-1)
+        output.PointData.append(y1, "scalars")
+        output.PointData.SetActiveScalars("scalars")
 
         return 1
 
@@ -141,7 +109,7 @@ class PythonDnsDataReader(VTKPythonAlgorithmBase):
     def SetThickness(self, x):
         self.Modified()
 
-    def _read_data (self, fname):
+    def _read_data (self, fname, only_tags=False):
         import struct
         import sys
         import os
@@ -161,14 +129,21 @@ class PythonDnsDataReader(VTKPythonAlgorithmBase):
            #dummy = f1.read(4)
 
            #raw_field = f1.read(4+bSize*5+4)[4:-4]
-            V = np.fromfile(f1, dtype=np.double, count=(3*count)).reshape((3,nz,ny,nx))
-            s = np.fromfile(f1, dtype=np.double, count=(2*count)).reshape((2,nz,ny,nx))
+            if only_tags:
+                u = None
+                v = None
+                w = None
+                Y0 = None
+                Y1 = None
+            else:
+                V = np.fromfile(f1, dtype=np.double, count=(3*count)).reshape((3,nz,ny,nx))
+                s = np.fromfile(f1, dtype=np.double, count=(2*count)).reshape((2,nz,ny,nx))
 
-            u = V[0]
-            v = V[1]
-            w = V[2]
+                u = V[0]
+                v = V[1]
+                w = V[2]
 
-            Y0 = s[0]
-            Y1 = s[1]
+                Y0 = s[0]
+                Y1 = s[1]
 
         return t, nx, ny, nz, u, v, w, Y0, Y1