diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/data_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/data_reader.py
index 1568b838..c931fdfc 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/data_reader.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/data_reader.py
@@ -138,7 +138,7 @@ def read_openPMD_params(self, iteration, extract_parameters=True):
                     self.series, iteration, extract_parameters)
 
     def read_field_cartesian( self, iteration, field, coord, axis_labels,
-                          slice_relative_position, slice_across ):
+                          slice_relative_position, slice_across, units ):
         """
         Extract a given field from an openPMD file in the openPMD format,
         when the geometry is cartesian (1d, 2d or 3d).
@@ -172,6 +172,9 @@ def read_field_cartesian( self, iteration, field, coord, axis_labels,
            0 : middle of the simulation box
            1 : upper edge of the simulation box
 
+        units: string
+           Type of units to be used for data reading.
+
         Returns
         -------
         A tuple with
@@ -183,14 +186,14 @@ def read_field_cartesian( self, iteration, field, coord, axis_labels,
             filename = self.iteration_to_file[iteration]
             return h5py_reader.read_field_cartesian(
                 filename, iteration, field, coord, axis_labels,
-                slice_relative_position, slice_across )
+                slice_relative_position, slice_across, units )
         elif self.backend == 'openpmd-api':
             return io_reader.read_field_cartesian(
                 self.series, iteration, field, coord, axis_labels,
-                slice_relative_position, slice_across )
+                slice_relative_position, slice_across, units )
 
     def read_field_circ( self, iteration, field, coord, slice_relative_position,
-                        slice_across, m=0, theta=0. ):
+                        slice_across, units, m=0, theta=0. ):
         """
         Extract a given field from an openPMD file in the openPMD format,
         when the geometry is thetaMode
@@ -229,6 +232,9 @@ def read_field_circ( self, iteration, field, coord, slice_relative_position,
            0 : middle of the simulation box
            1 : upper edge of the simulation box
 
+        units: string
+           Type of units to be used for data reading.
+
         Returns
         -------
         A tuple with
@@ -241,13 +247,14 @@ def read_field_circ( self, iteration, field, coord, slice_relative_position,
             filename = self.iteration_to_file[iteration]
             return h5py_reader.read_field_circ(
                 filename, iteration, field, coord, slice_relative_position,
-                slice_across, m, theta )
+                slice_across, units, m, theta )
         elif self.backend == 'openpmd-api':
             return io_reader.read_field_circ(
                 self.series, iteration, field, coord, slice_relative_position,
-                slice_across, m, theta )
+                slice_across, units, m, theta )
 
-    def read_species_data( self, iteration, species, record_comp, extensions):
+    def read_species_data( self, iteration, species, record_comp,
+                           extensions, units):
         """
         Extract a given species' record_comp
 
@@ -265,14 +272,19 @@ def read_species_data( self, iteration, species, record_comp, extensions):
 
         extensions: list of strings
             The extensions that the current OpenPMDTimeSeries complies with
+
+        units: string
+            Type of units to be used for data reading.
         """
         if self.backend == 'h5py':
             filename = self.iteration_to_file[iteration]
             return h5py_reader.read_species_data(
-                    filename, iteration, species, record_comp, extensions )
+                    filename, iteration, species, record_comp,
+                    extensions, units)
         elif self.backend == 'openpmd-api':
             return io_reader.read_species_data(
-                    self.series, iteration, species, record_comp, extensions )
+                    self.series, iteration, species, record_comp,
+                    extensions, units)
 
     def get_grid_parameters(self, iteration, avail_fields, metadata ):
         """
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/field_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/field_reader.py
index 43a8b0f9..301a033f 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/field_reader.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/field_reader.py
@@ -16,7 +16,7 @@
 
 
 def read_field_cartesian( filename, iteration, field, coord, axis_labels,
-                          slice_relative_position, slice_across ):
+                          slice_relative_position, slice_across, units ):
     """
     Extract a given field from an HDF5 file in the openPMD format,
     when the geometry is cartesian (1d, 2d or 3d).
@@ -53,6 +53,9 @@ def read_field_cartesian( filename, iteration, field, coord, axis_labels,
        0 : middle of the simulation box
        1 : upper edge of the simulation box
 
+    units: string
+       Type of units to be used for data reading.
+
     Returns
     -------
     A tuple with
@@ -103,11 +106,11 @@ def read_field_cartesian( filename, iteration, field, coord, axis_labels,
 
         axes = { i: axis_labels[i] for i in range(len(axis_labels)) }
         # Extract data
-        F = get_data( dset, list_i_cell, list_slicing_index )
+        F = get_data( dset, units, list_i_cell, list_slicing_index )
         info = FieldMetaInformation( axes, shape, grid_spacing, global_offset,
                 group.attrs['gridUnitSI'], dset.attrs['position'] )
     else:
-        F = get_data( dset )
+        F = get_data( dset, units )
         axes = { i: axis_labels[i] for i in range(len(axis_labels)) }
         info = FieldMetaInformation( axes, F.shape,
             group.attrs['gridSpacing'], group.attrs['gridGlobalOffset'],
@@ -119,7 +122,8 @@ def read_field_cartesian( filename, iteration, field, coord, axis_labels,
 
 
 def read_field_circ( filename, iteration, field, coord,
-                     slice_relative_position, slice_across, m=0, theta=0. ):
+                     slice_relative_position, slice_across,
+                     units, m=0, theta=0. ):
     """
     Extract a given field from an HDF5 file in the openPMD format,
     when the geometry is thetaMode
@@ -159,6 +163,9 @@ def read_field_circ( filename, iteration, field, coord,
        0 : middle of the simulation box
        1 : upper edge of the simulation box
 
+    units: string
+       Type of units to be used for data reading.
+
     Returns
     -------
     A tuple with
@@ -188,7 +195,7 @@ def read_field_circ( filename, iteration, field, coord,
         # Get cylindrical info
         rmax = info.rmax
         inv_dr = 1./info.dr
-        Fcirc = get_data( dset )  # (Extracts all modes)
+        Fcirc = get_data( dset, units )  # (Extracts all modes)
         nr = Fcirc.shape[1]
         if m == 'all':
             modes = [ mode for mode in range(0, int(Nm / 2) + 1) ]
@@ -221,22 +228,22 @@ def read_field_circ( filename, iteration, field, coord,
             mult_above_axis = np.array( mult_above_axis )
             mult_below_axis = np.array( mult_below_axis )
             # - Sum the modes
-            F = get_data( dset )  # (Extracts all modes)
+            F = get_data( dset, units )  # (Extracts all modes)
             F_total[Nr:, :] = np.tensordot( mult_above_axis,
                                             F, axes=(0, 0) )[:, :]
             F_total[:Nr, :] = np.tensordot( mult_below_axis,
                                             F, axes=(0, 0) )[::-1, :]
         elif m == 0:
             # Extract mode 0
-            F = get_data( dset, 0, 0 )
+            F = get_data( dset, units, 0, 0 )
             F_total[Nr:, :] = F[:, :]
             F_total[:Nr, :] = F[::-1, :]
         else:
             # Extract higher mode
             cos = np.cos( m * theta )
             sin = np.sin( m * theta )
-            F_cos = get_data( dset, 2 * m - 1, 0 )
-            F_sin = get_data( dset, 2 * m, 0 )
+            F_cos = get_data( dset, units, 2 * m - 1, 0 )
+            F_sin = get_data( dset, units, 2 * m, 0 )
             F = cos * F_cos + sin * F_sin
             F_total[Nr:, :] = F[:, :]
             F_total[:Nr, :] = (-1) ** m * F[::-1, :]
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/particle_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/particle_reader.py
index 82fddd20..3728edfe 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/particle_reader.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/particle_reader.py
@@ -15,7 +15,8 @@
 from .utilities import get_data, join_infile_path
 
 
-def read_species_data(filename, iteration, species, record_comp, extensions):
+def read_species_data(filename, iteration, species, record_comp,
+                      extensions, units):
     """
     Extract a given species' record_comp
 
@@ -36,6 +37,10 @@ def read_species_data(filename, iteration, species, record_comp, extensions):
 
     extensions: list of strings
         The extensions that the current OpenPMDTimeSeries complies with
+
+    units: string
+        Type of units to be used for data reading. Will convert ux uy uz to
+        normalized units if 'SI_u'
     """
     # Open the HDF5 file
     dfile = h5py.File( filename, 'r' )
@@ -63,7 +68,8 @@ def read_species_data(filename, iteration, species, record_comp, extensions):
         output_type = np.uint64
     else:
         output_type = np.float64
-    data = get_data( species_grp[ opmd_record_comp ], output_type=output_type )
+    data = get_data( species_grp[ opmd_record_comp ], units,
+                     output_type=output_type)
 
     # For ED-PIC: if the data is weighted for a full macroparticle,
     # divide by the weight with the proper power
@@ -74,16 +80,17 @@ def read_species_data(filename, iteration, species, record_comp, extensions):
         macro_weighted = record_dset.attrs['macroWeighted']
         weighting_power = record_dset.attrs['weightingPower']
         if (macro_weighted == 1) and (weighting_power != 0):
-            w = get_data( species_grp[ 'weighting' ] )
+            w = get_data( species_grp[ 'weighting' ], units )
             data *= w ** (-weighting_power)
 
     # - Return positions, with an offset
     if record_comp in ['x', 'y', 'z']:
-        offset = get_data(species_grp['positionOffset/%s' % record_comp])
+        offset = get_data(species_grp['positionOffset/%s' % record_comp],
+                          units)
         data += offset
     # - Return momentum in normalized units
-    elif record_comp in ['ux', 'uy', 'uz' ]:
-        m = get_data(species_grp['mass'])
+    elif record_comp in ['ux', 'uy', 'uz' ] and units == 'SI_u':
+        m = get_data(species_grp['mass'], units)
         # Normalize only if the particle mass is non-zero
         if np.all( m != 0 ):
             norm_factor = 1. / (m * constants.c)
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/utilities.py b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/utilities.py
index 5713a173..c1b14706 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/utilities.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/h5py_reader/utilities.py
@@ -75,7 +75,7 @@ def is_scalar_record(record):
     return(scalar)
 
 
-def get_data(dset, i_slice=None, pos_slice=None, output_type=np.float64):
+def get_data(dset, units, i_slice=None, pos_slice=None, output_type=np.float64):
     """
     Extract the data from a (possibly constant) dataset
     Slice the data according to the parameters i_slice and pos_slice
@@ -85,6 +85,10 @@ def get_data(dset, i_slice=None, pos_slice=None, output_type=np.float64):
     dset: an h5py.Dataset or h5py.Group (when constant)
         The object from which the data is extracted
 
+    units: string
+        Type of units to be used for data reading. Won't multiply data
+        with unit_SI when 'raw'
+
     pos_slice: int or list of int, optional
         Slice direction(s).
         When None, no slicing is performed
@@ -140,9 +144,8 @@ def get_data(dset, i_slice=None, pos_slice=None, output_type=np.float64):
         data = data.astype( output_type )
     # Scale by the conversion factor
     if output_type in [ np.float64, np.float32, np.float16 ]:
-        if dset.attrs['unitSI'] != 1.0:
-            data *= dset.attrs['unitSI']
-
+        if dset.attrs['unitSI'] != 1.0 and not units == 'raw':
+           data *= dset.attrs['unitSI']
     return(data)
 
 
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/field_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/field_reader.py
index 3b9447d4..c4b0302e 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/field_reader.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/field_reader.py
@@ -15,7 +15,8 @@
 
 
 def read_field_cartesian( series, iteration, field_name, component_name,
-                          axis_labels, slice_relative_position, slice_across ):
+                          axis_labels, slice_relative_position, slice_across,
+                          units):
     """
     Extract a given field from a file in the openPMD format,
     when the geometry is cartesian (1d, 2d or 3d).
@@ -52,6 +53,9 @@ def read_field_cartesian( series, iteration, field_name, component_name,
        0 : middle of the simulation box
        1 : upper edge of the simulation box
 
+    units: string
+       Type of units to be used for data reading.
+
     Returns
     -------
     A tuple with
@@ -106,21 +110,23 @@ def read_field_cartesian( series, iteration, field_name, component_name,
 
         axes = { i: axis_labels[i] for i in range(len(axis_labels)) }
         # Extract data
-        F = get_data( series, component, list_i_cell, list_slicing_index )
+        F = get_data( series, component, units,
+                      list_i_cell, list_slicing_index )
         info = FieldMetaInformation( axes, shape, grid_spacing, global_offset,
                 grid_unit_SI, grid_position )
     else:
-        F = get_data( series, component )
+        F = get_data( series, component, units)
         axes = { i: axis_labels[i] for i in range(len(axis_labels)) }
         info = FieldMetaInformation( axes, F.shape,
             grid_spacing, global_offset,
             grid_unit_SI, grid_position )
 
-    return F, info 
+    return F, info
 
 
 def read_field_circ( series, iteration, field_name, component_name,
-                     slice_relative_position, slice_across, m=0, theta=0. ):
+                     slice_relative_position, slice_across,
+                     units, m=0, theta=0. ):
     """
     Extract a given field from a file in the openPMD format,
     when the geometry is thetaMode
@@ -160,6 +166,9 @@ def read_field_circ( series, iteration, field_name, component_name,
        0 : middle of the simulation box
        1 : upper edge of the simulation box
 
+    units: string
+       Type of units to be used for data reading.
+
     Returns
     -------
     A tuple with
@@ -191,7 +200,7 @@ def read_field_circ( series, iteration, field_name, component_name,
         # Get cylindrical info
         rmax = info.rmax
         inv_dr = 1./info.dr
-        Fcirc = get_data( series, component )  # (Extracts all modes)
+        Fcirc = get_data( series, component, units)  # (Extracts all modes)
         nr = Fcirc.shape[1]
         if m == 'all':
             modes = [ mode for mode in range(0, int(Nm / 2) + 1) ]
@@ -224,22 +233,22 @@ def read_field_circ( series, iteration, field_name, component_name,
             mult_above_axis = np.array( mult_above_axis )
             mult_below_axis = np.array( mult_below_axis )
             # - Sum the modes
-            F = get_data( series, component )  # (Extracts all modes)
+            F = get_data( series, component, units )  # (Extracts all modes)
             F_total[Nr:, :] = np.tensordot( mult_above_axis,
                                             F, axes=(0, 0) )[:, :]
             F_total[:Nr, :] = np.tensordot( mult_below_axis,
                                             F, axes=(0, 0) )[::-1, :]
         elif m == 0:
             # Extract mode 0
-            F = get_data( series, component, 0, 0 )
+            F = get_data( series, component, units, 0, 0 )
             F_total[Nr:, :] = F[:, :]
             F_total[:Nr, :] = F[::-1, :]
         else:
             # Extract higher mode
             cos = np.cos( m * theta )
             sin = np.sin( m * theta )
-            F_cos = get_data( series, component, 2 * m - 1, 0 )
-            F_sin = get_data( series, component, 2 * m, 0 )
+            F_cos = get_data( series, component, units, 2 * m - 1, 0 )
+            F_sin = get_data( series, component, units, 2 * m, 0 )
             F = cos * F_cos + sin * F_sin
             F_total[Nr:, :] = F[:, :]
             F_total[:Nr, :] = (-1) ** m * F[::-1, :]
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/particle_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/particle_reader.py
index 0f14e866..fdd77ea8 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/particle_reader.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/particle_reader.py
@@ -15,7 +15,7 @@
 
 
 def read_species_data(series, iteration, species_name, component_name,
-                      extensions):
+                      extensions, units):
     """
     Extract a given species' record_comp
 
@@ -36,6 +36,10 @@ def read_species_data(series, iteration, species_name, component_name,
 
     extensions: list of strings
         The extensions that the current OpenPMDTimeSeries complies with
+
+    units: string
+        Type of units to be used for data reading. Will convert ux uy uz to
+        narmalized units if 'SI_u'
     """
     it = series.iterations[iteration]
 
@@ -47,7 +51,7 @@ def read_species_data(series, iteration, species_name, component_name,
                         'uy': ['momentum', 'y'],
                         'uz': ['momentum', 'z'],
                         'w': ['weighting', None]}
-    
+
     if component_name in dict_record_comp:
         ompd_record_name, ompd_record_comp_name = \
             dict_record_comp[component_name]
@@ -70,7 +74,7 @@ def read_species_data(series, iteration, species_name, component_name,
         output_type = np.uint64
     else:
         output_type = np.float64
-    data = get_data( series, component, output_type=output_type )
+    data = get_data( series, component, units, output_type=output_type)
 
     # For ED-PIC: if the data is weighted for a full macroparticle,
     # divide by the weight with the proper power
@@ -80,17 +84,18 @@ def read_species_data(series, iteration, species_name, component_name,
         weighting_power = record.get_attribute('weightingPower')
         if (macro_weighted == 1) and (weighting_power != 0):
             w_component = next(species['weighting'].items())[1]
-            w = get_data( w_component )
+            w = get_data(series, w_component, units)
             data *= w ** (-weighting_power)
 
     # - Return positions, with an offset
     if component_name in ['x', 'y', 'z']:
-        offset = get_data(series, species['positionOffset'][component_name])
+        offset = get_data(series, species['positionOffset'][component_name],
+                          units)
         data += offset
     # - Return momentum in normalized units
-    elif component_name in ['ux', 'uy', 'uz' ]:
+    elif component_name in ['ux', 'uy', 'uz' ] and units == 'SI_u':
         mass_component = next(species['mass'].items())[1]
-        m = get_data(series, mass_component)
+        m = get_data(series, mass_component, units)
         # Normalize only if the particle mass is non-zero
         if np.all( m != 0 ):
             norm_factor = 1. / (m * constants.c)
diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/utilities.py b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/utilities.py
index bb933287..6a3e97ad 100644
--- a/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/utilities.py
+++ b/openpmd_viewer/openpmd_timeseries/data_reader/io_reader/utilities.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 
-def get_data(series, record_component, i_slice=None, pos_slice=None,
+def get_data(series, record_component, units, i_slice=None, pos_slice=None,
              output_type=np.float64):
     """
     Extract the data from a (possibly constant) dataset
@@ -24,6 +24,10 @@ def get_data(series, record_component, i_slice=None, pos_slice=None,
 
     record_component: an openPMD.Record_Component
 
+    units: string
+        Type of units to be used for data reading. Won't multiply data
+        with unit_SI when 'raw'
+
     pos_slice: int or list of int, optional
         Slice direction(s).
         When None, no slicing is performed
@@ -70,9 +74,8 @@ def get_data(series, record_component, i_slice=None, pos_slice=None,
         data = data.astype( output_type )
     # Scale by the conversion factor
     if output_type in [ np.float64, np.float32, np.float16 ]:
-        if record_component.unit_SI != 1.0:
+        if record_component.unit_SI != 1.0 and not units == 'raw':
             data *= record_component.unit_SI
-
     return data
 
 
diff --git a/openpmd_viewer/openpmd_timeseries/main.py b/openpmd_viewer/openpmd_timeseries/main.py
index 3de3d9a3..185d5010 100644
--- a/openpmd_viewer/openpmd_timeseries/main.py
+++ b/openpmd_viewer/openpmd_timeseries/main.py
@@ -34,7 +34,8 @@ class OpenPMDTimeSeries(InteractiveViewer):
     - slider
     """
 
-    def __init__(self, path_to_dir, check_all_files=True, backend=None):
+    def __init__(self, path_to_dir, check_all_files=True, backend=None,
+                 units='SI_u'):
         """
         Initialize an openPMD time series
 
@@ -56,6 +57,11 @@ def __init__(self, path_to_dir, check_all_files=True, backend=None):
             Backend to be used for data reading. Can be `openpmd-api`
             or `h5py`. If not provided will use `openpmd-api` if available
             and `h5py` otherwise.
+
+        units: string
+            Type of units to be used for data reading. 'SI' for SI units,
+            'SI_u' for SI units but with momentum normalized by m_species * c,
+            'raw' to ignore the units_SI atribute for all components.
         """
         # Check backend
         if backend is None:
@@ -66,6 +72,12 @@ def __init__(self, path_to_dir, check_all_files=True, backend=None):
                     .format(backend, available_backends) )
         self.backend = backend
 
+        # Check if units is vailid
+        if units not in ['SI_u', 'SI', 'raw']:
+            raise OpenPMDException("Invalid value for argument units. "
+                "Must be 'SI_u', 'SI', or 'raw'")
+        self.units = units
+
         # Initialize data reader
         self.data_reader = DataReader(backend)
 
@@ -272,14 +284,15 @@ def get_particle(self, var_list=None, species=None, t=None, iteration=None,
         data_list = []
         for quantity in var_list:
             data_list.append( self.data_reader.read_species_data(
-                iteration, species, quantity, self.extensions))
+                iteration, species, quantity, self.extensions, self.units))
         # Apply selection if needed
         if isinstance( select, dict ):
             data_list = apply_selection( iteration, self.data_reader,
-                data_list, select, species, self.extensions)
+                data_list, select, species, self.extensions, self.units)
         elif isinstance( select, ParticleTracker ):
             data_list = select.extract_tracked_particles( iteration,
-                self.data_reader, data_list, species, self.extensions )
+                self.data_reader, data_list, species,
+                self.extensions, self.units)
 
         # Plotting
         if plot and len(var_list) in [1, 2]:
@@ -287,13 +300,14 @@ def get_particle(self, var_list=None, species=None, t=None, iteration=None,
             # Extract the weights, if they are available
             if 'w' in self.avail_record_components[species]:
                 w = self.data_reader.read_species_data(
-                    iteration, species, 'w', self.extensions)
+                    iteration, species, 'w', self.extensions, self.units)
                 if isinstance( select, dict ):
                     w, = apply_selection( iteration, self.data_reader,
-                        [w], select, species, self.extensions)
+                        [w], select, species, self.extensions, self.units)
                 elif isinstance( select, ParticleTracker ):
                     w, = select.extract_tracked_particles( iteration,
-                        self.data_reader, [w], species, self.extensions )
+                        self.data_reader, [w], species,
+                        self.extensions, self.units)
             # Otherwise consider that all particles have a weight of 1
             else:
                 w = np.ones_like(data_list[0])
@@ -502,7 +516,7 @@ def get_field(self, field=None, coord=None, t=None, iteration=None,
         if geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]:
             F, info = self.data_reader.read_field_cartesian(
                 iteration, field, coord, axis_labels,
-                slice_relative_position, slice_across)
+                slice_relative_position, slice_across, self.units)
         # - For thetaMode
         elif geometry == "thetaMode":
             if (coord in ['x', 'y']) and \
@@ -510,16 +524,16 @@ def get_field(self, field=None, coord=None, t=None, iteration=None,
                 # For Cartesian components, combine r and t components
                 Fr, info = self.data_reader.read_field_circ(
                     iteration, field, 'r', slice_relative_position,
-                    slice_across, m, theta)
+                    slice_across, self.units, m, theta)
                 Ft, info = self.data_reader.read_field_circ(
                     iteration, field, 't', slice_relative_position,
-                    slice_across, m, theta)
+                    slice_across, self.units, m, theta)
                 F = combine_cylindrical_components(Fr, Ft, theta, coord, info)
             else:
                 # For cylindrical or scalar components, no special treatment
                 F, info = self.data_reader.read_field_circ(iteration,
                     field, coord, slice_relative_position,
-                    slice_across, m, theta)
+                    slice_across, self.units, m, theta)
 
         # Plot the resulting field
         # Deactivate plotting when there is no slice selection
diff --git a/openpmd_viewer/openpmd_timeseries/particle_tracker.py b/openpmd_viewer/openpmd_timeseries/particle_tracker.py
index b51e6d98..74ea5a19 100644
--- a/openpmd_viewer/openpmd_timeseries/particle_tracker.py
+++ b/openpmd_viewer/openpmd_timeseries/particle_tracker.py
@@ -106,7 +106,7 @@ def __init__(self, ts, species=None, t=None,
 
 
     def extract_tracked_particles( self, iteration, data_reader, data_list,
-                                    species, extensions ):
+                                    species, extensions, units):
         """
         Select the elements of each particle quantities in data_list,
         so as to only return those that correspond to the tracked particles
@@ -129,6 +129,9 @@ def extract_tracked_particles( self, iteration, data_reader, data_list,
         extensions: list of strings
             The extensions that the current OpenPMDTimeSeries complies with
 
+        units: string
+            Type of units to be used for data reading.
+
         Returns
         -------
         A list of 1darrays that correspond to data_list, but where only the
@@ -137,7 +140,8 @@ def extract_tracked_particles( self, iteration, data_reader, data_list,
         initialization)
         """
         # Extract the particle id, and get the extraction indices
-        pid = data_reader.read_species_data(iteration, species, 'id', extensions)
+        pid = data_reader.read_species_data(iteration, species, 'id',
+                                            extensions, units)
         selected_indices = self.get_extraction_indices( pid )
 
         # For each particle quantity, select only the tracked particles
diff --git a/openpmd_viewer/openpmd_timeseries/utilities.py b/openpmd_viewer/openpmd_timeseries/utilities.py
index 6286ec8f..1e5240e6 100644
--- a/openpmd_viewer/openpmd_timeseries/utilities.py
+++ b/openpmd_viewer/openpmd_timeseries/utilities.py
@@ -52,7 +52,7 @@ def sanitize_slicing(slice_across, slice_relative_position):
     return copy.copy(slice_across), copy.copy(slice_relative_position)
 
 def apply_selection(iteration, data_reader, data_list,
-                    select, species, extensions):
+                    select, species, extensions, units):
     """
     Select the elements of each particle quantities in data_list,
     based on the selection rules in `select`
@@ -94,7 +94,7 @@ def apply_selection(iteration, data_reader, data_list,
     # Loop through the selection rules, and aggregate results in select_array
     for quantity in select.keys():
         q = data_reader.read_species_data(
-            iteration, species, quantity, extensions)
+            iteration, species, quantity, extensions, units)
         # Check lower bound
         if select[quantity][0] is not None:
             select_array = np.logical_and(