diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst
index d1862881588..47a7a9cb76f 100644
--- a/doc/sphinx/analysis.rst
+++ b/doc/sphinx/analysis.rst
@@ -345,8 +345,10 @@ Both versions are equivalent in the :math:`N\rightarrow \infty` limit but give n
 Observables and correlators
 ---------------------------
 
-Analysis in the core is a new concept introduced in since version 3.1.
-It was motivated by the fact, that sometimes it is desirable that the
+Observables extract properties of the particles and the LB fluid and
+return either the raw data or a statistic derived from them. The
+Observable framework is progressively replacing the Analysis framework.
+This is motivated by the fact, that sometimes it is desirable that the
 analysis functions do more than just return a value to the scripting
 interface. For some observables it is desirable to be sampled every few
 integration steps. In addition, it should be possible to pass the
@@ -380,12 +382,13 @@ not all observables carry out their calculations in parallel.
 Instead, the entire particle configuration is collected on the head node, and the calculations are carried out there.
 This is only performance-relevant if the number of processor cores is large and/or interactions are calculated very frequently.
 
-.. _Creating an observable:
+.. _Using observables:
 
-Creating an observable
-~~~~~~~~~~~~~~~~~~~~~~
+Using observables
+~~~~~~~~~~~~~~~~~
 
-The observables are represented as Python classes derived from :class:`espressomd.observables.Observable`. They are contained in
+The observables are represented as Python classes derived from
+:class:`espressomd.observables.Observable`. They are contained in
 the ``espressomd.observables`` module. An observable is instantiated as
 follows
 
@@ -397,17 +400,70 @@ follows
 Here, the keyword argument ``ids`` specifies the ids of the particles,
 which the observable should take into account.
 
-The current value of an observable can be obtained using its calculate()-method::
+The current value of an observable can be obtained using its
+:meth:`~espressomd.observables.Observable.calculate` method::
 
     print(part_pos.calculate())
 
+Profile observables have additional methods
+:meth:`~espressomd.observables.ProfileObservable.bin_centers` and
+:meth:`~espressomd.observables.ProfileObservable.bin_edges` to facilitate
+plotting of histogram slices with functions that require either bin centers
+or bin edges for the axes. Example::
+
+    import matplotlib.pyplot as plt
+    import numpy as np
+    import espressomd
+    import espressomd.observables
+
+    system = espressomd.System(box_l=[10.0, 10.0, 10.0])
+    system.part.add(id=0, pos=[4.0, 3.0, 6.0])
+    system.part.add(id=1, pos=[7.0, 3.0, 6.0])
+
+    # histogram in Cartesian coordinates
+    density_profile = espressomd.observables.DensityProfile(
+        ids=[0, 1],
+        n_x_bins=8, min_x=1.0, max_x=9.0,
+        n_y_bins=8, min_y=1.0, max_y=9.0,
+        n_z_bins=4, min_z=4.0, max_z=8.0)
+    obs_data = density_profile.calculate()
+    obs_bins = density_profile.bin_centers()
+
+    # 1D slice: requires bin centers
+    plt.plot(obs_bins[:, 2, 2, 0], obs_data[:, 2, 2])
+    plt.show()
+
+    # 2D slice: requires extent
+    plt.imshow(obs_data[:, :, 2].T, origin='lower',
+               extent=[density_profile.min_x, density_profile.max_x,
+                       density_profile.min_y, density_profile.max_y])
+    plt.show()
+
+    # histogram in cylindrical coordinates
+    density_profile = espressomd.observables.CylindricalDensityProfile(
+        ids=[0, 1], center=[5.0, 5.0, 0.0], axis=[0, 0, 1],
+        n_r_bins=8, min_r=0.0, max_r=4.0,
+        n_phi_bins=16, min_phi=-np.pi, max_phi=np.pi,
+        n_z_bins=4, min_z=4.0, max_z=8.0)
+    obs_data = density_profile.calculate()
+    obs_bins = density_profile.bin_edges()
+
+    # 2D slice: requires bin edges
+    fig = plt.figure()
+    ax = fig.add_subplot(111, polar='True')
+    r = obs_bins[:, 0, 0, 0]
+    phi = obs_bins[0, :, 0, 1]
+    ax.pcolormesh(phi, r, obs_data[:, :, 2])
+    plt.show()
+
+
 .. _Available observables:
 
 Available observables
 ~~~~~~~~~~~~~~~~~~~~~
 
-The following list contains some of the available observables. You can find documentation for
-all available observables in :mod:`espressomd.observables`.
+The following list contains some of the available observables. You can find
+documentation for all available observables in :mod:`espressomd.observables`.
 
 - Observables working on a given set of particles:
 
@@ -468,7 +524,6 @@ all available observables in :mod:`espressomd.observables`.
 
    - :class:`~espressomd.observables.CylindricalLBVelocityProfileAtParticlePositions`
 
-
 - System-wide observables
 
    - :class:`~espressomd.observables.StressTensor`: Total stress tensor (see :ref:`stress tensor`)
@@ -476,7 +531,6 @@ all available observables in :mod:`espressomd.observables`.
    - :class:`~espressomd.observables.DPDStress`
 
 
-
 .. _Correlations:
 
 Correlations
diff --git a/samples/lb_profile.py b/samples/lb_profile.py
index dba07d8a31b..9d599af91fa 100644
--- a/samples/lb_profile.py
+++ b/samples/lb_profile.py
@@ -44,7 +44,7 @@
 system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23)
 fluid_obs = espressomd.observables.CylindricalLBVelocityProfile(
     center=[5.0, 5.0, 0.0],
-    axis='z',
+    axis=[0, 0, 1],
     n_r_bins=100,
     n_phi_bins=1,
     n_z_bins=1,
@@ -54,9 +54,7 @@
     max_phi=np.pi,
     min_z=0.0,
     max_z=10.0,
-    sampling_delta_x=0.05,
-    sampling_delta_y=0.05,
-    sampling_delta_z=1.0)
+    sampling_density=0.1)
 cylinder_shape = espressomd.shapes.Cylinder(
     center=[5.0, 5.0, 5.0],
     axis=[0, 0, 1],
diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp
index db64b084534..241ebc7bbde 100644
--- a/src/core/observables/CylindricalLBProfileObservable.hpp
+++ b/src/core/observables/CylindricalLBProfileObservable.hpp
@@ -19,8 +19,7 @@
 #ifndef OBSERVABLES_CYLINDRICALLBPROFILEOBSERVABLE_HPP
 #define OBSERVABLES_CYLINDRICALLBPROFILEOBSERVABLE_HPP
 
-#include "CylindricalProfile.hpp"
-#include "Observable.hpp"
+#include "CylindricalProfileObservable.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/math/coordinate_transformation.hpp>
@@ -31,16 +30,16 @@ using Utils::Vector3d;
 
 namespace Observables {
 
-class CylindricalLBProfileObservable : public Observable,
-                                       public CylindricalProfile {
+class CylindricalLBProfileObservable : public CylindricalProfileObservable {
 public:
   CylindricalLBProfileObservable(Vector3d const &center, Vector3d const &axis,
                                  int n_r_bins, int n_phi_bins, int n_z_bins,
                                  double min_r, double min_phi, double min_z,
                                  double max_r, double max_phi, double max_z,
                                  double sampling_density)
-      : CylindricalProfile(center, axis, min_r, max_r, min_phi, max_phi, min_z,
-                           max_z, n_r_bins, n_phi_bins, n_z_bins),
+      : CylindricalProfileObservable(center, axis, min_r, max_r, min_phi,
+                                     max_phi, min_z, max_z, n_r_bins,
+                                     n_phi_bins, n_z_bins),
         sampling_density(sampling_density) {
     calculate_sampling_positions();
   }
diff --git a/src/core/observables/CylindricalPidProfileObservable.hpp b/src/core/observables/CylindricalPidProfileObservable.hpp
index 60932aef263..ffc576a54bb 100644
--- a/src/core/observables/CylindricalPidProfileObservable.hpp
+++ b/src/core/observables/CylindricalPidProfileObservable.hpp
@@ -19,15 +19,13 @@
 #ifndef OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP
 #define OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP
 
-#include <cmath>
-
-#include "CylindricalProfile.hpp"
+#include "CylindricalProfileObservable.hpp"
 #include "PidObservable.hpp"
 
 namespace Observables {
 
 class CylindricalPidProfileObservable : public PidObservable,
-                                        public CylindricalProfile {
+                                        public CylindricalProfileObservable {
 public:
   CylindricalPidProfileObservable(std::vector<int> const &ids,
                                   Utils::Vector3d const &center,
@@ -36,8 +34,9 @@ class CylindricalPidProfileObservable : public PidObservable,
                                   double min_phi, double min_z, double max_r,
                                   double max_phi, double max_z)
       : PidObservable(ids),
-        CylindricalProfile(center, axis, min_r, max_r, min_phi, max_phi, min_z,
-                           max_z, n_r_bins, n_phi_bins, n_z_bins) {}
+        CylindricalProfileObservable(center, axis, min_r, max_r, min_phi,
+                                     max_phi, min_z, max_z, n_r_bins,
+                                     n_phi_bins, n_z_bins) {}
 };
 
 } // Namespace Observables
diff --git a/src/core/observables/CylindricalProfile.hpp b/src/core/observables/CylindricalProfileObservable.hpp
similarity index 51%
rename from src/core/observables/CylindricalProfile.hpp
rename to src/core/observables/CylindricalProfileObservable.hpp
index a87e5e2506f..e443d9796ec 100644
--- a/src/core/observables/CylindricalProfile.hpp
+++ b/src/core/observables/CylindricalProfileObservable.hpp
@@ -19,17 +19,24 @@
 #ifndef OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP
 #define OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP
 
-#include <utils/Vector.hpp>
+#include <boost/range/algorithm.hpp>
+#include <vector>
+
+#include "Observable.hpp"
 
-#include <string>
+#include <utils/Vector.hpp>
+#include <utils/math/make_lin_space.hpp>
 
 namespace Observables {
-class CylindricalProfile {
+
+/** Cylindrical profile observable */
+class CylindricalProfileObservable : virtual public Observable {
 public:
-  CylindricalProfile(Utils::Vector3d const &center, Utils::Vector3d const &axis,
-                     double min_r, double max_r, double min_phi, double max_phi,
-                     double min_z, double max_z, int n_r_bins, int n_phi_bins,
-                     int n_z_bins)
+  CylindricalProfileObservable(Utils::Vector3d const &center,
+                               Utils::Vector3d const &axis, double min_r,
+                               double max_r, double min_phi, double max_phi,
+                               double min_z, double max_z, int n_r_bins,
+                               int n_phi_bins, int n_z_bins)
       : center(center), axis(axis), min_r(min_r), max_r(max_r),
         min_phi(min_phi), max_phi(max_phi), min_z(min_z), max_z(max_z),
         n_r_bins(static_cast<size_t>(n_r_bins)),
@@ -37,11 +44,30 @@ class CylindricalProfile {
         n_z_bins(static_cast<size_t>(n_z_bins)){};
   Utils::Vector3d center;
   Utils::Vector3d axis;
+  // Range of the profile edges.
   double min_r, max_r;
   double min_phi, max_phi;
   double min_z, max_z;
   // Number of bins for each coordinate.
   size_t n_r_bins, n_phi_bins, n_z_bins;
+
+  std::vector<size_t> shape() const override {
+    return {n_r_bins, n_phi_bins, n_z_bins};
+  }
+
+  /** Calculate the bin edges for each dimension */
+  std::array<std::vector<double>, 3> edges() {
+    std::array<std::vector<double>, 3> profile_edges = {
+        {std::vector<double>(n_r_bins + 1), std::vector<double>(n_phi_bins + 1),
+         std::vector<double>(n_z_bins + 1)}};
+    boost::copy(Utils::make_lin_space(min_r, max_r, n_r_bins + 1),
+                profile_edges[0].begin());
+    boost::copy(Utils::make_lin_space(min_phi, max_phi, n_phi_bins + 1),
+                profile_edges[1].begin());
+    boost::copy(Utils::make_lin_space(min_z, max_z, n_z_bins + 1),
+                profile_edges[2].begin());
+    return profile_edges;
+  }
 };
 
 } // Namespace Observables
diff --git a/src/core/observables/ProfileObservable.hpp b/src/core/observables/ProfileObservable.hpp
index f9626c836fc..81d1b459c9d 100644
--- a/src/core/observables/ProfileObservable.hpp
+++ b/src/core/observables/ProfileObservable.hpp
@@ -19,11 +19,16 @@
 #ifndef OBSERVABLES_PROFILEOBSERVABLE_HPP
 #define OBSERVABLES_PROFILEOBSERVABLE_HPP
 
+#include <boost/range/algorithm.hpp>
+#include <vector>
+
+#include <utils/math/make_lin_space.hpp>
+
 #include "Observable.hpp"
 
 namespace Observables {
 
-// Observable which acts on a given list of particle ids
+/** Cartesian profile observable */
 class ProfileObservable : virtual public Observable {
 public:
   ProfileObservable(double min_x, double max_x, double min_y, double max_y,
@@ -33,13 +38,30 @@ class ProfileObservable : virtual public Observable {
         max_z(max_z), n_x_bins(static_cast<size_t>(n_x_bins)),
         n_y_bins(static_cast<size_t>(n_y_bins)),
         n_z_bins(static_cast<size_t>(n_z_bins)) {}
+  // Range of the profile edges.
   double min_x, max_x;
   double min_y, max_y;
   double min_z, max_z;
+  // Number of bins for each coordinate.
   size_t n_x_bins, n_y_bins, n_z_bins;
+
   std::vector<size_t> shape() const override {
     return {n_x_bins, n_y_bins, n_z_bins};
   }
+
+  /** Calculate the bin edges for each dimension */
+  std::array<std::vector<double>, 3> edges() {
+    std::array<std::vector<double>, 3> profile_edges = {
+        {std::vector<double>(n_x_bins + 1), std::vector<double>(n_y_bins + 1),
+         std::vector<double>(n_z_bins + 1)}};
+    boost::copy(Utils::make_lin_space(min_x, max_x, n_x_bins + 1),
+                profile_edges[0].begin());
+    boost::copy(Utils::make_lin_space(min_y, max_y, n_y_bins + 1),
+                profile_edges[1].begin());
+    boost::copy(Utils::make_lin_space(min_z, max_z, n_z_bins + 1),
+                profile_edges[2].begin());
+    return profile_edges;
+  }
 };
 
 } // Namespace Observables
diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py
index 9dbc23b460d..8ff0cc98384 100644
--- a/src/python/espressomd/observables.py
+++ b/src/python/espressomd/observables.py
@@ -14,12 +14,21 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+import itertools
 import numpy as np
 from .script_interface import ScriptInterfaceHelper, script_interface_register
 
 
 @script_interface_register
 class Observable(ScriptInterfaceHelper):
+    """
+    Base class for all observables.
+
+    Methods
+    -------
+    shape()
+        Return the shape of the observable.
+    """
     _so_name = "Observables::Observable"
     _so_bind_methods = ("shape",)
     _so_creation_policy = "LOCAL"
@@ -28,6 +37,38 @@ def calculate(self):
         return np.array(self.call_method("calculate")).reshape(self.shape())
 
 
+class ProfileObservable(Observable):
+    """
+    Base class for histogram-based observables.
+    """
+
+    def bin_edges(self):
+        """
+        Returns
+        -------
+        :obj:`ndarray` of :obj:`float`
+            Positions between the bins. If the histogram has dimensions
+            ``(M,N,O)``, the bin edges have dimensions ``(M+1,N+1,O+1,3)``.
+        """
+        edges = self.call_method("edges")
+        shape = list(map(len, edges)) + [len(edges)]
+        return np.array(list(itertools.product(*edges))).reshape(shape)
+
+    def bin_centers(self):
+        """
+        Returns
+        -------
+        :obj:`ndarray` of :obj:`float`
+            Positions of the bins centers. If the histogram has dimensions
+            ``(M,N,O)``, the bin centers have dimensions ``(M,N,O,3)``.
+        """
+        edges = self.call_method("edges")
+        for i, edge in enumerate(edges):
+            edges[i] = np.array(edge[:-1]) + (edge[1] - edge[0]) / 2
+        shape = list(map(len, edges)) + [len(edges)]
+        return np.array(list(itertools.product(*edges))).reshape(shape)
+
+
 @script_interface_register
 class ComPosition(Observable):
 
@@ -93,7 +134,7 @@ class Current(Observable):
 
 
 @script_interface_register
-class DensityProfile(Observable):
+class DensityProfile(ProfileObservable):
 
     """Calculates the particle density profile for particles with given ids.
 
@@ -149,7 +190,7 @@ class DipoleMoment(Observable):
 
 
 @script_interface_register
-class FluxDensityProfile(Observable):
+class FluxDensityProfile(ProfileObservable):
 
     """Calculates the particle flux density for particles with given ids.
 
@@ -187,7 +228,7 @@ class FluxDensityProfile(Observable):
 
 
 @script_interface_register
-class ForceDensityProfile(Observable):
+class ForceDensityProfile(ProfileObservable):
 
     """Calculates the force density profile for particles with given ids.
 
@@ -225,7 +266,7 @@ class ForceDensityProfile(Observable):
 
 
 @script_interface_register
-class LBVelocityProfile(Observable):
+class LBVelocityProfile(ProfileObservable):
 
     """Calculates the LB fluid velocity profile.
 
@@ -589,7 +630,7 @@ class DPDStress(Observable):
 
 
 @script_interface_register
-class CylindricalDensityProfile(Observable):
+class CylindricalDensityProfile(ProfileObservable):
 
     """Calculates the particle density in cylindrical coordinates.
 
@@ -629,7 +670,7 @@ class CylindricalDensityProfile(Observable):
 
 
 @script_interface_register
-class CylindricalFluxDensityProfile(Observable):
+class CylindricalFluxDensityProfile(ProfileObservable):
 
     """Calculates the particle flux density in cylindrical coordinates.
 
@@ -671,7 +712,7 @@ class CylindricalFluxDensityProfile(Observable):
 
 
 @script_interface_register
-class CylindricalLBFluxDensityProfileAtParticlePositions(Observable):
+class CylindricalLBFluxDensityProfileAtParticlePositions(ProfileObservable):
 
     """Calculates the LB fluid flux density at the particle positions in
     cylindrical coordinates.
@@ -714,7 +755,7 @@ class CylindricalLBFluxDensityProfileAtParticlePositions(Observable):
 
 
 @script_interface_register
-class CylindricalLBVelocityProfileAtParticlePositions(Observable):
+class CylindricalLBVelocityProfileAtParticlePositions(ProfileObservable):
 
     """Calculates the LB fluid velocity at the particle positions in
     cylindrical coordinates.
@@ -757,7 +798,7 @@ class CylindricalLBVelocityProfileAtParticlePositions(Observable):
 
 
 @script_interface_register
-class CylindricalVelocityProfile(Observable):
+class CylindricalVelocityProfile(ProfileObservable):
 
     """Calculates the particle velocity profile in cylindrical coordinates.
 
@@ -799,7 +840,7 @@ class CylindricalVelocityProfile(Observable):
 
 
 @script_interface_register
-class CylindricalLBVelocityProfile(Observable):
+class CylindricalLBVelocityProfile(ProfileObservable):
 
     """Calculates the LB fluid velocity profile in cylindrical coordinates.
 
diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp
index 95847f79147..6bceb55e14e 100644
--- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp
+++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp
@@ -24,6 +24,8 @@
 
 #include "script_interface/auto_parameters/AutoParameters.hpp"
 
+#include <boost/range/algorithm.hpp>
+#include <iterator>
 #include <memory>
 
 #include "Observable.hpp"
@@ -37,10 +39,14 @@ template <typename CoreCylLBObs>
 class CylindricalLBProfileObservable
     : public AutoParameters<CylindricalLBProfileObservable<CoreCylLBObs>,
                             Observable> {
+  using Base =
+      AutoParameters<CylindricalLBProfileObservable<CoreCylLBObs>, Observable>;
+
 public:
   static_assert(std::is_base_of<::Observables::CylindricalLBProfileObservable,
                                 CoreCylLBObs>::value,
                 "");
+  using Base::Base;
   CylindricalLBProfileObservable() {
     this->add_parameters({
         {"center",
@@ -133,14 +139,13 @@ class CylindricalLBProfileObservable
 
   Variant call_method(std::string const &method,
                       VariantMap const &parameters) override {
-    if (method == "calculate") {
-      return cylindrical_profile_observable()->operator()();
-    }
-    if (method == "shape") {
-      auto const shape = cylindrical_profile_observable()->shape();
-      return std::vector<int>{shape.begin(), shape.end()};
+    if (method == "edges") {
+      std::vector<Variant> variant_edges;
+      boost::copy(cylindrical_profile_observable()->edges(),
+                  std::back_inserter(variant_edges));
+      return variant_edges;
     }
-    return {};
+    return Base::call_method(method, parameters);
   }
 
   std::shared_ptr<::Observables::Observable> observable() const override {
diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp
index 4ae01fbef83..4ffd84aa5cb 100644
--- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp
+++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp
@@ -22,6 +22,8 @@
 #ifndef SCRIPT_INTERFACE_OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP
 #define SCRIPT_INTERFACE_OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP
 
+#include <boost/range/algorithm.hpp>
+#include <iterator>
 #include <memory>
 
 #include "script_interface/ScriptInterface.hpp"
@@ -37,10 +39,14 @@ template <typename CoreObs>
 class CylindricalPidProfileObservable
     : public AutoParameters<CylindricalPidProfileObservable<CoreObs>,
                             Observable> {
+  using Base =
+      AutoParameters<CylindricalPidProfileObservable<CoreObs>, Observable>;
+
 public:
   static_assert(std::is_base_of<::Observables::CylindricalPidProfileObservable,
                                 CoreObs>::value,
                 "");
+  using Base::Base;
   CylindricalPidProfileObservable() {
     this->add_parameters({
         {"ids",
@@ -131,6 +137,17 @@ class CylindricalPidProfileObservable
             "max_z");
   }
 
+  Variant call_method(std::string const &method,
+                      VariantMap const &parameters) override {
+    if (method == "edges") {
+      std::vector<Variant> variant_edges;
+      boost::copy(cylindrical_pid_profile_observable()->edges(),
+                  std::back_inserter(variant_edges));
+      return variant_edges;
+    }
+    return Base::call_method(method, parameters);
+  }
+
   std::shared_ptr<::Observables::Observable> observable() const override {
     return m_observable;
   }
diff --git a/src/script_interface/observables/LBProfileObservable.hpp b/src/script_interface/observables/LBProfileObservable.hpp
index 626d1deb25b..c616f71a204 100644
--- a/src/script_interface/observables/LBProfileObservable.hpp
+++ b/src/script_interface/observables/LBProfileObservable.hpp
@@ -24,6 +24,8 @@
 
 #include "script_interface/auto_parameters/AutoParameters.hpp"
 
+#include <boost/range/algorithm.hpp>
+#include <iterator>
 #include <memory>
 
 #include "Observable.hpp"
@@ -35,10 +37,13 @@ namespace Observables {
 template <typename CoreLBObs>
 class LBProfileObservable
     : public AutoParameters<LBProfileObservable<CoreLBObs>, Observable> {
+  using Base = AutoParameters<LBProfileObservable<CoreLBObs>, Observable>;
+
 public:
   static_assert(
       std::is_base_of<::Observables::LBProfileObservable, CoreLBObs>::value,
       "");
+  using Base::Base;
   LBProfileObservable() {
     this->add_parameters(
         {{"n_x_bins",
@@ -145,14 +150,13 @@ class LBProfileObservable
 
   Variant call_method(std::string const &method,
                       VariantMap const &parameters) override {
-    if (method == "calculate") {
-      return profile_observable()->operator()();
-    }
-    if (method == "shape") {
-      auto const shape = profile_observable()->shape();
-      return std::vector<int>{shape.begin(), shape.end()};
+    if (method == "edges") {
+      std::vector<Variant> variant_edges;
+      boost::copy(profile_observable()->edges(),
+                  std::back_inserter(variant_edges));
+      return variant_edges;
     }
-    return {};
+    return Base::call_method(method, parameters);
   }
 
   std::shared_ptr<::Observables::Observable> observable() const override {
diff --git a/src/script_interface/observables/PidProfileObservable.hpp b/src/script_interface/observables/PidProfileObservable.hpp
index 499bf72c722..1a817d8cb76 100644
--- a/src/script_interface/observables/PidProfileObservable.hpp
+++ b/src/script_interface/observables/PidProfileObservable.hpp
@@ -23,6 +23,9 @@
 #define SCRIPT_INTERFACE_OBSERVABLES_PIDPROFILEOBSERVABLE_HPP
 
 #include "script_interface/auto_parameters/AutoParameters.hpp"
+
+#include <boost/range/algorithm.hpp>
+#include <iterator>
 #include <memory>
 
 #include "core/observables/DensityProfile.hpp"
@@ -36,7 +39,10 @@ namespace Observables {
 template <typename CoreObs>
 class PidProfileObservable
     : public AutoParameters<PidProfileObservable<CoreObs>, Observable> {
+  using Base = AutoParameters<PidProfileObservable<CoreObs>, Observable>;
+
 public:
+  using Base::Base;
   PidProfileObservable() {
     this->add_parameters(
         {{"ids",
@@ -108,6 +114,17 @@ class PidProfileObservable
             "min_z", "max_x", "max_y", "max_z");
   }
 
+  Variant call_method(std::string const &method,
+                      VariantMap const &parameters) override {
+    if (method == "edges") {
+      std::vector<Variant> variant_edges;
+      boost::copy(pid_profile_observable()->edges(),
+                  std::back_inserter(variant_edges));
+      return variant_edges;
+    }
+    return Base::call_method(method, parameters);
+  }
+
   std::shared_ptr<::Observables::PidProfileObservable>
   pid_profile_observable() const {
     return m_observable;
diff --git a/src/script_interface/observables/ProfileObservable.hpp b/src/script_interface/observables/ProfileObservable.hpp
index dc464385003..6916b6f5a43 100644
--- a/src/script_interface/observables/ProfileObservable.hpp
+++ b/src/script_interface/observables/ProfileObservable.hpp
@@ -34,7 +34,10 @@ namespace Observables {
 template <typename CoreObs>
 class ProfileObservable
     : public AutoParameters<ProfileObservable<CoreObs>, Observable> {
+  using Base = AutoParameters<ProfileObservable<CoreObs>, Observable>;
+
 public:
+  using Base::Base;
   ProfileObservable() {
     this->add_parameters(
         {{"n_x_bins",
@@ -95,6 +98,17 @@ class ProfileObservable
 
   void construct(VariantMap const &params) override {}
 
+  Variant call_method(std::string const &method,
+                      VariantMap const &parameters) override {
+    if (method == "edges") {
+      std::vector<Variant> variant_edges;
+      boost::copy(profile_observable()->edges(),
+                  std::back_inserter(variant_edges));
+      return variant_edges;
+    }
+    return Base::call_method(method, parameters);
+  }
+
   std::shared_ptr<::Observables::ProfileObservable> profile_observable() const {
     return m_observable;
   }
diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py
index 3915c07f3dc..3be0875408e 100644
--- a/testsuite/python/observable_cylindrical.py
+++ b/testsuite/python/observable_cylindrical.py
@@ -14,7 +14,6 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-import sys
 import numpy as np
 import unittest as ut
 import espressomd
@@ -28,9 +27,8 @@ class TestCylindricalObservable(ut.TestCase):
     Testcase for the cylindrical observables.
 
     """
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system = espressomd.System(box_l=[15.0, 15.0, 15.0])
     system.time_step = 0.01
-    system.box_l = [15.0, 15.0, 15.0]
     system.cell_system.skin = 0.4
 
     params = {
@@ -104,15 +102,9 @@ def set_particles(self):
 
     def calculate_numpy_histogram(self):
         pol_positions = self.pol_coords()
-        np_hist, _ = np.histogramdd(
-            pol_positions,
-            bins=(self.params['n_r_bins'],
-                  self.params['n_phi_bins'],
-                  self.params['n_z_bins']),
-            range=[(self.params['min_r'], self.params['max_r']),
-                   (self.params['min_phi'], self.params['max_phi']),
-                   (self.params['min_z'], self.params['max_z'])])
-        return np_hist
+        np_hist, np_edges = tests_common.get_histogram(
+            pol_positions, self.params, 'cylindrical')
+        return np_hist, np_edges
 
     def normalize_with_bin_volume(self, histogram):
         bin_volume = tests_common.get_cylindrical_bin_volume(
@@ -141,9 +133,12 @@ def density_profile_test(self):
             local_params['axis'] = [0.0, 0.0, 1.0]
         obs = espressomd.observables.CylindricalDensityProfile(**local_params)
         core_hist = obs.calculate()
-        np_hist = self.calculate_numpy_histogram()
+        core_edges = obs.call_method("edges")
+        np_hist, np_edges = self.calculate_numpy_histogram()
         np_hist = self.normalize_with_bin_volume(np_hist)
         np.testing.assert_array_almost_equal(np_hist, core_hist)
+        for i in range(3):
+            np.testing.assert_array_almost_equal(np_edges[i], core_edges[i])
         self.assertEqual(np.prod(obs.shape()), len(np_hist.flatten()))
 
     def velocity_profile_test(self):
@@ -161,7 +156,7 @@ def velocity_profile_test(self):
         core_hist_v_r = core_hist[:, :, :, 0]
         core_hist_v_phi = core_hist[:, :, :, 1]
         core_hist_v_z = core_hist[:, :, :, 2]
-        np_hist = self.calculate_numpy_histogram()
+        np_hist, _ = self.calculate_numpy_histogram()
         for x in np.nditer(np_hist, op_flags=['readwrite']):
             if x[...] > 0.0:
                 x[...] /= x[...]
@@ -187,7 +182,7 @@ def flux_density_profile_test(self):
         core_hist_v_r = core_hist[:, :, :, 0]
         core_hist_v_phi = core_hist[:, :, :, 1]
         core_hist_v_z = core_hist[:, :, :, 2]
-        np_hist = self.calculate_numpy_histogram()
+        np_hist, _ = self.calculate_numpy_histogram()
         np_hist = self.normalize_with_bin_volume(np_hist)
         np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r)
         np.testing.assert_array_almost_equal(
@@ -215,8 +210,4 @@ def test_hist_z(self):
 
 
 if __name__ == "__main__":
-    suite = ut.TestSuite()
-    suite.addTests(ut.TestLoader().loadTestsFromTestCase(
-        TestCylindricalObservable))
-    result = ut.TextTestRunner(verbosity=4).run(suite)
-    sys.exit(not result.wasSuccessful())
+    ut.main()
diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py
index e955808f548..55b9e8caed7 100644
--- a/testsuite/python/observable_cylindricalLB.py
+++ b/testsuite/python/observable_cylindricalLB.py
@@ -147,23 +147,17 @@ def LB_fluxdensity_profile_test(self):
         core_hist_v_r = core_hist[:, :, :, 0]
         core_hist_v_phi = core_hist[:, :, :, 1]
         core_hist_v_z = core_hist[:, :, :, 2]
+        core_edges = p.call_method("edges")
         self.pol_positions = self.pol_coords()
-        np_hist, _ = np.histogramdd(
-            self.pol_positions,
-            bins=(self.params['n_r_bins'],
-                  self.params['n_phi_bins'],
-                  self.params['n_z_bins']),
-            range=[(self.params['min_r'],
-                    self.params['max_r']),
-                   (self.params['min_phi'],
-                    self.params['max_phi']),
-                   (self.params['min_z'],
-                    self.params['max_z'])])
+        np_hist, np_edges = tests_common.get_histogram(
+            self.pol_positions, self.params, 'cylindrical')
         np_hist = self.normalize_with_bin_volume(np_hist)
         np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r)
         np.testing.assert_array_almost_equal(
             np_hist * self.v_phi, core_hist_v_phi)
         np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z)
+        for i in range(3):
+            np.testing.assert_array_almost_equal(np_edges[i], core_edges[i])
         self.assertEqual(np.prod(p.shape()), len(np_hist.flatten()) * 3)
 
     def LB_velocity_profile_at_particle_positions_test(self):
diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py
index d28548f2a03..4661bf65c54 100644
--- a/testsuite/python/observable_profile.py
+++ b/testsuite/python/observable_profile.py
@@ -14,17 +14,16 @@
 #
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
-import sys
 import unittest as ut
 import unittest_decorators as utx
 import numpy as np
 import espressomd
 import espressomd.observables
+import tests_common
 
 
 class ProfileObservablesTest(ut.TestCase):
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
-    system.box_l = [10.0, 10.0, 10.0]
+    system = espressomd.System(box_l=[10.0, 10.0, 10.0])
     system.cell_system.skin = 0.1
     system.time_step = 0.01
     system.part.add(id=0, pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0])
@@ -44,9 +43,26 @@ class ProfileObservablesTest(ut.TestCase):
     def test_density_profile(self):
         density_profile = espressomd.observables.DensityProfile(**self.kwargs)
         obs_data = density_profile.calculate()
-        self.assertEqual(obs_data[0, 0, 1], 2.0 / self.bin_volume)
+        obs_edges = density_profile.call_method("edges")
+        obs_bin_edges = density_profile.bin_edges()
+        obs_bin_centers = density_profile.bin_centers()
+        np_hist, np_edges = tests_common.get_histogram(
+            np.copy(self.system.part[:].pos), self.kwargs, 'cartesian',
+            normed=True)
+        np_hist *= len(self.system.part)
+        np.testing.assert_array_almost_equal(obs_data, np_hist)
+        for i in range(3):
+            np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i])
         self.assertEqual(np.prod(obs_data.shape),
                          self.kwargs['n_x_bins'] * self.kwargs['n_y_bins'] * self.kwargs['n_z_bins'])
+        np.testing.assert_array_almost_equal(
+            obs_bin_edges[0, 0, 0],
+            [self.kwargs['min_x'], self.kwargs['min_y'], self.kwargs['min_z']])
+        np.testing.assert_array_almost_equal(
+            obs_bin_edges[-1, -1, -1],
+            [self.kwargs['max_x'], self.kwargs['max_y'], self.kwargs['max_z']])
+        np.testing.assert_array_almost_equal(obs_bin_centers[0, 0, 0],
+                                             [2.5, 2.5, 2.5])
 
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
     def test_force_density_profile(self):
@@ -71,8 +87,4 @@ def test_flux_density_profile(self):
 
 
 if __name__ == '__main__':
-    suite = ut.TestSuite()
-    suite.addTests(ut.TestLoader().loadTestsFromTestCase(
-        ProfileObservablesTest))
-    result = ut.TextTestRunner(verbosity=4).run(suite)
-    sys.exit(not result.wasSuccessful())
+    ut.main()
diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py
index 875b985e81d..fcfffb374eb 100644
--- a/testsuite/python/observable_profileLB.py
+++ b/testsuite/python/observable_profileLB.py
@@ -18,6 +18,7 @@
 import unittest_decorators as utx
 import numpy as np
 import copy
+import tests_common
 
 import espressomd
 import espressomd.lb
@@ -78,6 +79,12 @@ def test_velocity_profile(self):
         obs = espressomd.observables.LBVelocityProfile(
             **LB_VELOCITY_PROFILE_PARAMS)
         obs_data = obs.calculate()
+        obs_edges = obs.call_method("edges")
+        _, np_edges = tests_common.get_histogram(
+            np.zeros([1, 3]), LB_VELOCITY_PROFILE_PARAMS, 'cartesian',
+            normed=True)
+        for i in range(3):
+            np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i])
         for x in range(obs_data.shape[0]):
             for y in range(obs_data.shape[1]):
                 for z in range(obs_data.shape[2]):
diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py
index ce015a733f7..e685ec59220 100644
--- a/testsuite/python/tests_common.py
+++ b/testsuite/python/tests_common.py
@@ -275,6 +275,46 @@ def get_cylindrical_bin_volume(
             phi_bin_size / (2.0 * np.pi) * z_bin_size
     return bin_volume
 
+
+def get_histogram(pos, obs_params, coord_system, **kwargs):
+    """
+    Helper function for ``np.histogramdd()`` and observables.
+
+    Parameters
+    ----------
+    pos : (N, 3) array_like of :obj:`float`
+        Particle positions.
+    obs_params : :obj:`dict`
+        Parameters of the observable.
+    coord_system : :obj:`str`, \{'cartesian', 'cylindrical'\}
+        Coordinate system.
+    \*\*kwargs :
+        Optional parameters to ``np.histogramdd()``.
+
+    Returns
+    -------
+    array_like
+        Bins and bin edges.
+
+    """
+    if coord_system == 'cartesian':
+        bins = (obs_params['n_x_bins'],
+                obs_params['n_y_bins'],
+                obs_params['n_z_bins'])
+        extent = [(obs_params['min_x'], obs_params['max_x']),
+                  (obs_params['min_y'], obs_params['max_y']),
+                  (obs_params['min_z'], obs_params['max_z'])]
+    elif coord_system == 'cylindrical':
+        bins = (obs_params['n_r_bins'],
+                obs_params['n_phi_bins'],
+                obs_params['n_z_bins'])
+        extent = [(obs_params['min_r'], obs_params['max_r']),
+                  (obs_params['min_phi'], obs_params['max_phi']),
+                  (obs_params['min_z'], obs_params['max_z'])]
+    else:
+        raise ValueError("Unknown coord system '{}'".format(coord_system))
+    return np.histogramdd(pos, bins=bins, range=extent, **kwargs)
+
 #
 # Analytical Expressions for interactions
 #
diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt
index 09c4c87bb9a..2320b351b26 100644
--- a/testsuite/scripts/samples/CMakeLists.txt
+++ b/testsuite/scripts/samples/CMakeLists.txt
@@ -38,6 +38,7 @@ if(HDF5_FOUND)
 endif()
 sample_test(FILE test_lbf.py SUFFIX cpu)
 sample_test(FILE test_lbf.py SUFFIX gpu LABELS "gpu")
+sample_test(FILE test_lb_profile.py)
 sample_test(FILE test_lj_liquid_distribution.py)
 sample_test(FILE test_lj_liquid.py)
 sample_test(FILE test_lj_liquid_structurefactor.py)
diff --git a/testsuite/scripts/samples/test_lb_profile.py b/testsuite/scripts/samples/test_lb_profile.py
index 079d9339749..ecf91ade378 100644
--- a/testsuite/scripts/samples/test_lb_profile.py
+++ b/testsuite/scripts/samples/test_lb_profile.py
@@ -31,7 +31,7 @@ def test_fit(self):
         # compare the simulated data against the analytical solution
         sim = sample.lb_fluid_profile[:, 0, 0, 2]
         ana = sample.poiseuille_flow(sample.r, sample.r_max, 0.15)
-        self.assertLess(np.max(np.abs(sim - ana)), 0.12)
+        self.assertLess(np.max(np.abs(sim - ana)), 0.16)
 
 
 if __name__ == "__main__":