Repository: martinfleis/momepy
Pull request #78 (2019-09-16 14:42:21)

API: classes 
All functions measuring characters on gdf are now classes. This change enables better extensibility of each of them and provides a way of storing all used attributes and data.

codecov[bot]: # [Codecov](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=h1) Report
> Merging [#78](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=desc) into [master](https://codecov.io/gh/martinfleis/momepy/commit/a81ce021009dbb543aa226f5670f62c26c5bae7b?src=pr&el=desc) will **decrease** coverage by `0.06%`.
> The diff coverage is `96.54%`.

[![Impacted file tree graph](https://codecov.io/gh/martinfleis/momepy/pull/78/graphs/tree.svg?width=650&token=i04999UxkK&height=150&src=pr)](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master      #78      +/-   ##
==========================================
- Coverage   93.27%   93.21%   -0.07%     
==========================================
  Files           9        9              
  Lines        2336     2403      +67     
==========================================
+ Hits         2179     2240      +61     
- Misses        157      163       +6
```


| [Impacted Files](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [momepy/utils.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L3V0aWxzLnB5) | `81.52% <100%> (-0.1%)` | :arrow_down: |
| [momepy/diversity.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L2RpdmVyc2l0eS5weQ==) | `90.64% <90.32%> (+0.88%)` | :arrow_up: |
| [momepy/shape.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L3NoYXBlLnB5) | `95.21% <94.94%> (+0.03%)` | :arrow_up: |
| [momepy/intensity.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L2ludGVuc2l0eS5weQ==) | `97.51% <97.43%> (+0.35%)` | :arrow_up: |
| [momepy/dimension.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L2RpbWVuc2lvbi5weQ==) | `97.18% <98.26%> (-1.84%)` | :arrow_down: |
| [momepy/distribution.py](https://codecov.io/gh/martinfleis/momepy/pull/78/diff?src=pr&el=tree#diff-bW9tZXB5L2Rpc3RyaWJ1dGlvbi5weQ==) | `98.4% <98.37%> (+0.08%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=footer). Last update [a81ce02...4381eea](https://codecov.io/gh/martinfleis/momepy/pull/78?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/docs/api.rst b/docs/api.rst
index 45051ea..40fce3a 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -26,18 +26,18 @@ dimension
 .. autosummary::
    :toctree: generated/
 
-   area
-   average_character
-   courtyard_area
-   covered_area
-   floor_area
-   longest_axis_length
-   perimeter
-   segments_length
-   street_profile
-   volume
-   wall
-   weighted_character
+   Area
+   AverageCharacter
+   CourtyardArea
+   CoveredArea
+   FloorArea
+   LongestAxisLength
+   Perimeter
+   PerimeterWall
+   SegmentsLength
+   StreetProfile
+   Volume
+   WeightedCharacter
 
 shape
 -----
@@ -45,51 +45,51 @@ shape
 .. autosummary::
   :toctree: generated/
 
-  centroid_corners
-  circular_compactness
-  compactness_weighted_axis
-  convexeity
-  corners
-  courtyard_index
-  elongation
-  equivalent_rectangular_index
-  form_factor
-  fractal_dimension
-  linearity
-  rectangularity
-  shape_index
-  square_compactness
-  squareness
-  volume_facade_ratio
+  CentroidCorners
+  CircularCompactness
+  CompactnessWeightedAxis
+  Convexeity
+  Corners
+  CourtyardIndex
+  Elongation
+  EquivalentRectangularIndex
+  FormFactor
+  FractalDimension
+  Linearity
+  Rectangularity
+  ShapeIndex
+  SquareCompactness
+  Squareness
+  VolumeFacadeRatio
 
 spatial distribution
 --------------------
 .. autosummary::
   :toctree: generated/
 
-  alignment
-  building_adjacency
-  cell_alignment
-  mean_interbuilding_distance
-  neighbour_distance
-  neighbouring_street_orientation_deviation
-  neighbours
-  orientation
-  shared_walls_ratio
-  street_alignment
+  Alignment
+  BuildingAdjacency
+  CellAlignment
+  MeanInterbuildingDistance
+  NeighborDistance
+  NeighboringStreetOrientationDeviation
+  Neighbors
+  Orientation
+  SharedWallsRatio
+  StreetAlignment
 
 intensity
 ---------
 .. autosummary::
    :toctree: generated/
 
-   blocks_count
-   courtyards
-   density
-   elements_count
-   node_density
-   object_area_ratio
-   reached
+   AreaRatio
+   BlocksCount
+   Count
+   Courtyards
+   Density
+   NodeDensity
+   Reached
 
 graph
 ------------
@@ -118,10 +118,10 @@ diversity
 .. autosummary::
    :toctree: generated/
 
-   rng
-   simpson
-   theil
-   gini
+   Range
+   Simpson
+   Theil
+   Gini
 
 utilities
 ---------
diff --git a/docs/conf.py b/docs/conf.py
index dc91318..f30ced8 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -45,11 +45,14 @@ release = version
 # ones.
 extensions = [
     "sphinx.ext.autodoc",
-    "sphinx.ext.intersphinx",
     "sphinx.ext.autosummary",
     "sphinx.ext.viewcode",
-    "sphinx.ext.napoleon",
+    "sphinxcontrib.bibtex",
     "sphinx.ext.mathjax",
+    "sphinx.ext.doctest",
+    "sphinx.ext.intersphinx",
+    "numpydoc",
+    "matplotlib.sphinxext.plot_directive",
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -77,7 +80,7 @@ language = None
 exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
 
 # The name of the Pygments (syntax highlighting) style to use.
-pygments_style = None
+pygments_style = "sphinx"
 
 
 # -- Options for HTML output -------------------------------------------------
diff --git a/docs/environment.yml b/docs/environment.yml
index 6e91dbd..11cb520 100644
--- a/docs/environment.yml
+++ b/docs/environment.yml
@@ -15,3 +15,5 @@ dependencies:
     - osmnx
     - pip:
         - inequality
+        - sphinxcontrib.bibtex
+        - numpydoc
diff --git a/momepy/dimension.py b/momepy/dimension.py
index 4216f99..4e820e6 100644
--- a/momepy/dimension.py
+++ b/momepy/dimension.py
@@ -13,80 +13,76 @@ import numpy as np
 import scipy as sp
 
 
-def area(gdf):
+class Area:
     """
-    Calculates area of each object in given shapefile. It can be used for any
+    Calculates area of each object in given GeoDataFrame. It can be used for any
     suitable element (building footprint, plot, tessellation, block).
 
-    It is a simple wrapper for geopandas `gdf.geometry.area` for consistency of momepy.
+    It is a simple wrapper for geopandas `gdf.geometry.area` for the consistency of momepy.
 
     Parameters
     ----------
     gdf : GeoDataFrame
         GeoDataFrame containing objects to analyse
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    area : Series
+        Series containing resulting values
+
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     Examples
     --------
     >>> buildings = gpd.read_file(momepy.datasets.get_path('bubenec'), layer='buildings')
-    >>> buildings['area'] = momepy.area(buildings)
-    Calculating areas...
-    Areas calculated.
+    >>> buildings['area'] = momepy.Area(buildings).area
     >>> buildings.area[0]
     728.5574947044363
 
     """
 
-    print("Calculating areas...")
-
-    series = gdf.geometry.area
-
-    print("Areas calculated.")
-    return series
+    def __init__(self, gdf):
+        self.gdf = gdf
+        self.area = self.gdf.geometry.area
 
 
-def perimeter(gdf):
+class Perimeter:
     """
-    Calculates perimeter of each object in given shapefile. It can be used for any
+    Calculates perimeter of each object in given GeoDataFrame. It can be used for any
     suitable element (building footprint, plot, tessellation, block).
 
-    It is a simple wrapper for geopandas `gdf.geometry.length` for consistency of momepy.
+    It is a simple wrapper for geopandas `gdf.geometry.length` for the consistency of momepy.
 
     Parameters
     ----------
     gdf : GeoDataFrame
         GeoDataFrame containing objects to analyse
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    perimeter : Series
+        Series containing resulting values
+
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     Examples
     --------
     >>> buildings = gpd.read_file(momepy.datasets.get_path('bubenec'), layer='buildings')
-    >>> buildings['perimeter'] = momepy.perimeter(buildings)
-    Calculating perimeters...
-    Perimeters calculated.
+    >>> buildings['perimeter'] = momepy.Perimeter(buildings).perimeter
     >>> buildings.perimeter[0]
     137.18630991119903
     """
 
-    print("Calculating perimeters...")
-
-    series = gdf.geometry.length
+    def __init__(self, gdf):
+        self.gdf = gdf
+        self.perimeter = self.gdf.geometry.length
 
-    print("Perimeters calculated.")
-    return series
 
-
-def volume(gdf, heights, areas=None):
+class Volume:
     """
-    Calculates volume of each object in given shapefile based on its height and area.
+    Calculates volume of each object in given GeoDataFrame based on its height and area.
 
     .. math::
         area * height
@@ -101,50 +97,57 @@ def volume(gdf, heights, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    volume : Series
+        Series containing resulting values
+
+    gdf : GeoDataFrame
+        original GeoDataFrame
+
+    heights : Series
+        Series containing used heights values
+
+    areas : GeoDataFrame
+        Series containing used areas values
 
     Examples
     --------
-    >>> buildings['volume'] = momepy.volume(buildings, heights='height_col')
-    Calculating volumes...
-    Volumes calculated.
+    >>> buildings['volume'] = momepy.Volume(buildings, heights='height_col')
     >>> buildings.volume[0]
     7285.5749470443625
 
-    >>> buildings['volume'] = momepy.volume(buildings, heights='height_col', areas='area_col')
-    Calculating volumes...
-    Volumes calculated.
+    >>> buildings['volume'] = momepy.Volume(buildings, heights='height_col', areas='area_col').volume
     >>> buildings.volume[0]
     7285.5749470443625
     """
-    print("Calculating volumes...")
-    gdf = gdf.copy()
-    if not isinstance(heights, str):
-        gdf["mm_h"] = heights
-        heights = "mm_h"
 
-    if areas is not None:
-        if not isinstance(areas, str):
-            gdf["mm_a"] = areas
-            areas = "mm_a"
+    def __init__(self, gdf, heights, areas=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        if not isinstance(heights, str):
+            gdf["mm_h"] = heights
+            heights = "mm_h"
+        self.heights = gdf[heights]
+
+        if areas is not None:
+            if not isinstance(areas, str):
+                gdf["mm_a"] = areas
+                areas = "mm_a"
+            self.areas = gdf[areas]
+        else:
+            self.areas = gdf.geometry.area
         try:
-            series = gdf[areas] * gdf[heights]
+            self.volume = self.areas * self.heights
 
         except KeyError:
             raise KeyError(
                 "ERROR: Column not found. Define heights and areas or set areas to None."
             )
-    else:
-        series = gdf.geometry.area * gdf[heights]
 
-    print("Volumes calculated.")
-    return series
 
-
-def floor_area(gdf, heights, areas=None):
+class FloorArea:
     """
     Calculates floor area of each object based on height and area.
 
@@ -164,50 +167,56 @@ def floor_area(gdf, heights, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    fa : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    heights : Series
+        Series containing used heights values
+    areas : GeoDataFrame
+        Series containing used areas values
 
     Examples
     --------
-    >>> buildings['floor_area'] = momepy.floor_area(buildings, heights='height_col')
+    >>> buildings['floor_area'] = momepy.FloorArea(buildings, heights='height_col')
     Calculating floor areas...
     Floor areas calculated.
     >>> buildings.floor_area[0]
     2185.672484113309
 
-    >>> buildings['floor_area'] = momepy.floor_area(buildings, heights='height_col', areas='area_col')
-    Calculating floor areas...
-    Floor areas calculated.
+    >>> buildings['floor_area'] = momepy.FloorArea(buildings, heights='height_col', areas='area_col').fa
     >>> buildings.floor_area[0]
     2185.672484113309
     """
-    print("Calculating floor areas...")
-    gdf = gdf.copy()
-    if not isinstance(heights, str):
-        gdf["mm_h"] = heights
-        heights = "mm_h"
 
-    if areas is not None:
-        if not isinstance(areas, str):
-            gdf["mm_a"] = areas
-            areas = "mm_a"
+    def __init__(self, gdf, heights, areas=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        if not isinstance(heights, str):
+            gdf["mm_h"] = heights
+            heights = "mm_h"
+        self.heights = gdf[heights]
+
+        if areas is not None:
+            if not isinstance(areas, str):
+                gdf["mm_a"] = areas
+                areas = "mm_a"
+            self.areas = gdf[areas]
+        else:
+            self.areas = gdf.geometry.area
         try:
-            series = gdf[areas] * (gdf[heights] // 3)
+            self.fa = self.areas * (self.heights // 3)
 
         except KeyError:
             raise KeyError(
                 "ERROR: Column not found. Define heights and areas or set areas to None."
             )
-    else:
-        series = gdf.geometry.area * (gdf[heights] // 3)
 
-    print("Floor areas calculated.")
-    return series
 
-
-def courtyard_area(gdf, areas=None):
+class CourtyardArea:
     """
     Calculates area of holes within geometry - area of courtyards.
 
@@ -221,47 +230,42 @@ def courtyard_area(gdf, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ca : Series
+        Series containing resulting values
+
+    gdf : GeoDataFrame
+        original GeoDataFrame
+
+    areas : GeoDataFrame
+        Series containing used areas values
 
     Examples
     --------
-    >>> buildings['courtyard_area'] = momepy.courtyard_area(buildings)
-    Calculating courtyard areas...
-    Courtyard areas calculated.
+    >>> buildings['courtyard_area'] = momepy.CourtyardArea(buildings).ca
     >>> buildings.courtyard_area[80]
     353.33274206543274
     """
 
-    print("Calculating courtyard areas...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        if areas is None:
+            areas = gdf.geometry.area
 
-    if areas is not None:
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
-            lambda row: Polygon(row.geometry.exterior).area - row[areas], axis=1
-        )
+        self.areas = gdf[areas]
 
-    else:
-        series = gdf.apply(
-            lambda row: Polygon(row.geometry.exterior).area - row.geometry.area, axis=1
+        self.ca = gdf.apply(
+            lambda row: Polygon(row.geometry.exterior).area - row[areas], axis=1
         )
 
-    print("Courtyard areas calculated.")
-    return series
-
 
-# calculate the radius of circumcircle
-def _longest_axis(points):
-    circ = _make_circle(points)
-    return circ[2] * 2
-
-
-def longest_axis_length(gdf):
+class LongestAxisLength:
     """
     Calculates the length of the longest axis of object.
 
@@ -276,31 +280,35 @@ def longest_axis_length(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects to analyse
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    lal : Series
+        Series containing resulting values
+
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     Examples
     --------
-    >>> buildings['lal'] = momepy.longest_axis_length(buildings)
-    Calculating the longest axis...
-    The longest axis calculated.
+    >>> buildings['lal'] = momepy.LongestAxisLength(buildings).lal
     >>> buildings.lal[0]
     40.2655616057102
     """
 
-    print("Calculating the longest axis...")
-
-    series = gdf.apply(
-        lambda row: _longest_axis(row["geometry"].convex_hull.exterior.coords), axis=1
-    )
+    def __init__(self, gdf):
+        self.gdf = gdf
+        self.lal = gdf.apply(
+            lambda row: self._longest_axis(row.geometry.convex_hull.exterior.coords),
+            axis=1,
+        )
 
-    print("The longest axis calculated.")
-    return series
+    # calculate the radius of circumcircle
+    def _longest_axis(self, points):
+        circ = _make_circle(points)
+        return circ[2] * 2
 
 
-def average_character(gdf, values, spatial_weights, unique_id, rng=None, mode="mean"):
+class AverageCharacter:
     """
     Calculates the average of a character within k steps of morphological tessellation
 
@@ -326,10 +334,22 @@ def average_character(gdf, values, spatial_weights, unique_id, rng=None, mode="m
     mode : str (default 'mean')
         mode of average calculation. Can be set to `mean`, `median` or `mode`.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ac : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : GeoDataFrame
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    rng : tuple
+        range
+    mode : str
+        mode
 
     References
     ----------
@@ -339,52 +359,55 @@ def average_character(gdf, values, spatial_weights, unique_id, rng=None, mode="m
     Examples
     --------
     >>> sw = libpysal.weights.DistanceBand.from_dataframe(tessellation, threshold=100, silence_warnings=True, ids='uID')
-    >>> tessellation['mesh_100'] = momepy.average_character(tessellation, values='area', spatial_weights=sw, unique_id='uID')
-    Calculating average character value...
+    >>> tessellation['mesh_100'] = momepy.AverageCharacter(tessellation, values='area', spatial_weights=sw, unique_id='uID').ac
     100%|██████████| 144/144 [00:00<00:00, 1433.32it/s]
-    Average character value calculated.
     >>> tessellation.mesh_100[0]
     4823.1334436678835
     """
-    # define empty list for results
-    results_list = []
-
-    print("Calculating average character value...")
-    data = gdf.copy()
-
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
-
-    data = data.set_index(unique_id)
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        neighbours.append(index)
-
-        values_list = data.loc[neighbours][values]
+    def __init__(self, gdf, values, spatial_weights, unique_id, rng=None, mode="mean"):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.rng = rng
+        self.mode = mode
+
+        data = gdf.copy()
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
+
+        data = data.set_index(unique_id)
+
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
 
-        if rng:
-            from momepy import limit_range
+            values_list = data.loc[neighbours][values]
 
-            values_list = limit_range(values_list.tolist(), rng=rng)
-        if mode == "mean":
-            results_list.append(np.mean(values_list))
-        elif mode == "median":
-            results_list.append(np.median(values_list))
-        elif mode == "mode":
-            results_list.append(sp.stats.mode(values_list)[0][0])
-        else:
-            raise ValueError("{} is not supported as mode.".format(mode))
+            if rng:
+                from momepy import limit_range
 
-    series = pd.Series(results_list, index=gdf.index)
+                values_list = limit_range(values_list.tolist(), rng=rng)
+            if mode == "mean":
+                results_list.append(np.mean(values_list))
+            elif mode == "median":
+                results_list.append(np.median(values_list))
+            elif mode == "mode":
+                results_list.append(sp.stats.mode(values_list)[0][0])
+            else:
+                raise ValueError("{} is not supported as mode.".format(mode))
 
-    print("Average character value calculated.")
-    return series
+        self.ac = pd.Series(results_list, index=gdf.index)
 
 
-def street_profile(left, right, heights=None, distance=10, tick_length=50):
+class StreetProfile:
     """
     Calculates the street profile characters.
 
@@ -411,22 +434,30 @@ def street_profile(left, right, heights=None, distance=10, tick_length=50):
     tick_length : int (default 50)
         lenght of ticks
 
-    Returns
-    -------
-    street_profile : dictionary
-
-    'widths' : Series
-        Series containing street profile width values.
-    'width_deviations' : Series
-        Series containing street profile standard deviation values.
-    'openness' : Series
-        Series containing street profile openness values.
-    'heights' : Series, optional
+    Attributes
+    ----------
+    w : Series
+        Series containing street profile width values
+    wd : Series
+        Series containing street profile standard deviation values
+    o : Series
+        Series containing street profile openness values
+    h : Series
         Series containing street profile heights values. Returned only when heights is set.
-    'heights_deviations' : Series, optional
+    hd : Series
         Series containing street profile heights standard deviation values. Returned only when heights is set.
-    'profile' : Series, optional
+    p : Series
         Series containing street profile height/width ratio values. Returned only when heights is set.
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    distance : int
+        distance between perpendicular ticks
+    tick_length : int
+        lenght of ticks
+    heights : GeoDataFrame
+        Series containing used height values
 
     References
     ----------
@@ -439,27 +470,194 @@ def street_profile(left, right, heights=None, distance=10, tick_length=50):
 
     Examples
     --------
-    >>> street_profile = momepy.street_profile(streets_df, buildings_df, heights='height')
-    Calculating street profile...
+    >>> street_profile = momepy.StreetProfile(streets_df, buildings_df, heights='height')
     100%|██████████| 33/33 [00:02<00:00, 15.66it/s]
-    Street profile calculated.
-    >>> streets_df['width'] = street_profile['widths']
-    >>> streets_df['deviations'] = street_profile['width_deviations']
+    >>> streets_df['width'] = street_profile.w
+    >>> streets_df['deviations'] = street_profile.wd
     """
 
-    print("Calculating street profile...")
+    def __init__(self, left, right, heights=None, distance=10, tick_length=50):
+        self.left = left
+        self.right = right
+        self.distance = distance
+        self.tick_length = tick_length
+
+        if not isinstance(heights, str):
+            right = right.copy()
+            right["mm_h"] = heights
+            heights = "mm_h"
+
+        self.heights = right[heights]
+
+        sindex = right.sindex
+
+        results_list = []
+        deviations_list = []
+        heights_list = []
+        heights_deviations_list = []
+        openness_list = []
+
+        for idx, row in tqdm(left.iterrows(), total=left.shape[0]):
+            # if idx == 2:
+            #     ee
+            # list to hold all the point coords
+            list_points = []
+            # set the current distance to place the point
+            current_dist = distance
+            # make shapely MultiLineString object
+            shapely_line = row.geometry
+            # get the total length of the line
+            line_length = shapely_line.length
+            # append the starting coordinate to the list
+            list_points.append(Point(list(shapely_line.coords)[0]))
+            # https://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/
+            # while the current cumulative distance is less than the total length of the line
+            while current_dist < line_length:
+                # use interpolate and increase the current distance
+                list_points.append(shapely_line.interpolate(current_dist))
+                current_dist += distance
+            # append end coordinate to the list
+            list_points.append(Point(list(shapely_line.coords)[-1]))
+
+            ticks = []
+            for num, pt in enumerate(list_points, 1):
+                # start chainage 0
+                if num == 1:
+                    angle = self._getAngle(pt, list_points[num])
+                    line_end_1 = self._getPoint1(pt, angle, tick_length / 2)
+                    angle = self._getAngle(line_end_1, pt)
+                    line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
+                    tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
+                    tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
+                    ticks.append([tick1, tick2])
+
+                # everything in between
+                if num < len(list_points) - 1:
+                    angle = self._getAngle(pt, list_points[num])
+                    line_end_1 = self._getPoint1(
+                        list_points[num], angle, tick_length / 2
+                    )
+                    angle = self._getAngle(line_end_1, list_points[num])
+                    line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
+                    tick1 = LineString(
+                        [
+                            (line_end_1.x, line_end_1.y),
+                            (list_points[num].x, list_points[num].y),
+                        ]
+                    )
+                    tick2 = LineString(
+                        [
+                            (line_end_2.x, line_end_2.y),
+                            (list_points[num].x, list_points[num].y),
+                        ]
+                    )
+                    ticks.append([tick1, tick2])
+
+                # end chainage
+                if num == len(list_points):
+                    angle = self._getAngle(list_points[num - 2], pt)
+                    line_end_1 = self._getPoint1(pt, angle, tick_length / 2)
+                    angle = self._getAngle(line_end_1, pt)
+                    line_end_2 = self._getPoint2(line_end_1, angle, tick_length)
+                    tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
+                    tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
+                    ticks.append([tick1, tick2])
+            # widths = []
+            m_heights = []
+            lefts = []
+            rights = []
+            for duo in ticks:
+
+                for ix, tick in enumerate(duo):
+                    possible_intersections_index = list(
+                        sindex.intersection(tick.bounds)
+                    )
+                    possible_intersections = right.iloc[possible_intersections_index]
+                    real_intersections = possible_intersections.intersects(tick)
+                    get_height = right.loc[list(real_intersections.index)]
+                    possible_int = get_height.exterior.intersection(tick)
+
+                    if not possible_int.is_empty.all():
+                        true_int = []
+                        for one in list(possible_int.index):
+                            if possible_int[one].type == "Point":
+                                true_int.append(possible_int[one])
+                            elif possible_int[one].type == "MultiPoint":
+                                for p in possible_int[one]:
+                                    true_int.append(p)
+
+                        if len(true_int) > 1:
+                            distances = []
+                            ix = 0
+                            for p in true_int:
+                                distance = p.distance(Point(tick.coords[-1]))
+                                distances.append(distance)
+                                ix = ix + 1
+                            minimal = min(distances)
+                            if ix == 0:
+                                lefts.append(minimal)
+                            else:
+                                rights.append(minimal)
+                        else:
+                            if ix == 0:
+                                lefts.append(
+                                    true_int[0].distance(Point(tick.coords[-1]))
+                                )
+                            else:
+                                rights.append(
+                                    true_int[0].distance(Point(tick.coords[-1]))
+                                )
+                        if heights is not None:
+                            indices = {}
+                            for idx, row in get_height.iterrows():
+                                dist = row.geometry.distance(Point(tick.coords[-1]))
+                                indices[idx] = dist
+                            minim = min(indices, key=indices.get)
+                            m_heights.append(right.loc[minim][heights])
+
+            openness = (len(lefts) + len(rights)) / len(ticks * 2)
+            openness_list.append(1 - openness)
+            if rights and lefts:
+                results_list.append(2 * np.mean(lefts + rights))
+                deviations_list.append(np.std(lefts + rights))
+            elif not lefts and rights:
+                results_list.append(2 * np.mean([np.mean(rights), tick_length / 2]))
+                deviations_list.append(np.std(rights))
+            elif not rights and lefts:
+                results_list.append(2 * np.mean([np.mean(lefts), tick_length / 2]))
+                deviations_list.append(np.std(lefts))
+            else:
+                results_list.append(tick_length)
+                deviations_list.append(0)
+
+            if heights is not None:
+                if m_heights:
+                    heights_list.append(np.mean(m_heights))
+                    heights_deviations_list.append(np.std(m_heights))
+                else:
+                    heights_list.append(0)
+                    heights_deviations_list.append(0)
+
+        self.w = pd.Series(results_list, index=left.index)
+        self.wd = pd.Series(deviations_list, index=left.index)
+        self.o = pd.Series(openness_list, index=left.index)
+
+        if heights is not None:
+            self.h = pd.Series(heights_list, index=left.index)
+            self.hd = pd.Series(heights_deviations_list, index=left.index)
+            self.p = self.h / self.w
 
     # http://wikicode.wikidot.com/get-angle-of-line-between-two-points
     # https://glenbambrick.com/tag/perpendicular/
     # angle between two points
-    def _getAngle(pt1, pt2):
+    def _getAngle(self, pt1, pt2):
         x_diff = pt2.x - pt1.x
         y_diff = pt2.y - pt1.y
         return math.degrees(math.atan2(y_diff, x_diff))
 
     # start and end points of chainage tick
     # get the first end point of a tick
-    def _getPoint1(pt, bearing, dist):
+    def _getPoint1(self, pt, bearing, dist):
         angle = bearing + 90
         bearing = math.radians(angle)
         x = pt.x + dist * math.cos(bearing)
@@ -467,176 +665,14 @@ def street_profile(left, right, heights=None, distance=10, tick_length=50):
         return Point(x, y)
 
     # get the second end point of a tick
-    def _getPoint2(pt, bearing, dist):
+    def _getPoint2(self, pt, bearing, dist):
         bearing = math.radians(bearing)
         x = pt.x + dist * math.cos(bearing)
         y = pt.y + dist * math.sin(bearing)
         return Point(x, y)
 
-    sindex = right.sindex
-
-    results_list = []
-    deviations_list = []
-    heights_list = []
-    heights_deviations_list = []
-    openness_list = []
-
-    if heights is not None:
-        if not isinstance(heights, str):
-            right = right.copy()
-            right["mm_h"] = heights
-            heights = "mm_h"
-
-    for idx, row in tqdm(left.iterrows(), total=left.shape[0]):
-        # if idx == 2:
-        #     ee
-        # list to hold all the point coords
-        list_points = []
-        # set the current distance to place the point
-        current_dist = distance
-        # make shapely MultiLineString object
-        shapely_line = row.geometry
-        # get the total length of the line
-        line_length = shapely_line.length
-        # append the starting coordinate to the list
-        list_points.append(Point(list(shapely_line.coords)[0]))
-        # https://nathanw.net/2012/08/05/generating-chainage-distance-nodes-in-qgis/
-        # while the current cumulative distance is less than the total length of the line
-        while current_dist < line_length:
-            # use interpolate and increase the current distance
-            list_points.append(shapely_line.interpolate(current_dist))
-            current_dist += distance
-        # append end coordinate to the list
-        list_points.append(Point(list(shapely_line.coords)[-1]))
-
-        ticks = []
-        for num, pt in enumerate(list_points, 1):
-            # start chainage 0
-            if num == 1:
-                angle = _getAngle(pt, list_points[num])
-                line_end_1 = _getPoint1(pt, angle, tick_length / 2)
-                angle = _getAngle(line_end_1, pt)
-                line_end_2 = _getPoint2(line_end_1, angle, tick_length)
-                tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
-                tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
-                ticks.append([tick1, tick2])
-
-            # everything in between
-            if num < len(list_points) - 1:
-                angle = _getAngle(pt, list_points[num])
-                line_end_1 = _getPoint1(list_points[num], angle, tick_length / 2)
-                angle = _getAngle(line_end_1, list_points[num])
-                line_end_2 = _getPoint2(line_end_1, angle, tick_length)
-                tick1 = LineString(
-                    [
-                        (line_end_1.x, line_end_1.y),
-                        (list_points[num].x, list_points[num].y),
-                    ]
-                )
-                tick2 = LineString(
-                    [
-                        (line_end_2.x, line_end_2.y),
-                        (list_points[num].x, list_points[num].y),
-                    ]
-                )
-                ticks.append([tick1, tick2])
-
-            # end chainage
-            if num == len(list_points):
-                angle = _getAngle(list_points[num - 2], pt)
-                line_end_1 = _getPoint1(pt, angle, tick_length / 2)
-                angle = _getAngle(line_end_1, pt)
-                line_end_2 = _getPoint2(line_end_1, angle, tick_length)
-                tick1 = LineString([(line_end_1.x, line_end_1.y), (pt.x, pt.y)])
-                tick2 = LineString([(line_end_2.x, line_end_2.y), (pt.x, pt.y)])
-                ticks.append([tick1, tick2])
-        # widths = []
-        m_heights = []
-        lefts = []
-        rights = []
-        for duo in ticks:
-
-            for ix, tick in enumerate(duo):
-                possible_intersections_index = list(sindex.intersection(tick.bounds))
-                possible_intersections = right.iloc[possible_intersections_index]
-                real_intersections = possible_intersections.intersects(tick)
-                get_height = right.loc[list(real_intersections.index)]
-                possible_int = get_height.exterior.intersection(tick)
-
-                if not possible_int.is_empty.all():
-                    true_int = []
-                    for one in list(possible_int.index):
-                        if possible_int[one].type == "Point":
-                            true_int.append(possible_int[one])
-                        elif possible_int[one].type == "MultiPoint":
-                            for p in possible_int[one]:
-                                true_int.append(p)
-
-                    if len(true_int) > 1:
-                        distances = []
-                        ix = 0
-                        for p in true_int:
-                            distance = p.distance(Point(tick.coords[-1]))
-                            distances.append(distance)
-                            ix = ix + 1
-                        minimal = min(distances)
-                        if ix == 0:
-                            lefts.append(minimal)
-                        else:
-                            rights.append(minimal)
-                    else:
-                        if ix == 0:
-                            lefts.append(true_int[0].distance(Point(tick.coords[-1])))
-                        else:
-                            rights.append(true_int[0].distance(Point(tick.coords[-1])))
-                    if heights is not None:
-                        indices = {}
-                        for idx, row in get_height.iterrows():
-                            dist = row.geometry.distance(Point(tick.coords[-1]))
-                            indices[idx] = dist
-                        minim = min(indices, key=indices.get)
-                        m_heights.append(right.loc[minim][heights])
-
-        openness = (len(lefts) + len(rights)) / len(ticks * 2)
-        openness_list.append(1 - openness)
-        if rights and lefts:
-            results_list.append(2 * np.mean(lefts + rights))
-            deviations_list.append(np.std(lefts + rights))
-        elif not lefts and rights:
-            results_list.append(2 * np.mean([np.mean(rights), tick_length / 2]))
-            deviations_list.append(np.std(rights))
-        elif not rights and lefts:
-            results_list.append(2 * np.mean([np.mean(lefts), tick_length / 2]))
-            deviations_list.append(np.std(lefts))
-        else:
-            results_list.append(tick_length)
-            deviations_list.append(0)
-
-        if heights is not None:
-            if m_heights:
-                heights_list.append(np.mean(m_heights))
-                heights_deviations_list.append(np.std(m_heights))
-            else:
-                heights_list.append(0)
-                heights_deviations_list.append(0)
-
-    street_profile = {}
-    street_profile["widths"] = pd.Series(results_list, index=left.index)
-    street_profile["width_deviations"] = pd.Series(deviations_list, index=left.index)
-    street_profile["openness"] = pd.Series(openness_list, index=left.index)
-
-    if heights is not None:
-        street_profile["heights"] = pd.Series(heights_list, index=left.index)
-        street_profile["heights_deviations"] = pd.Series(
-            heights_deviations_list, index=left.index
-        )
-        street_profile["profile"] = street_profile["heights"] / street_profile["widths"]
 
-    print("Street profile calculated.")
-    return street_profile
-
-
-def weighted_character(gdf, values, spatial_weights, unique_id, areas=None):
+class WeightedCharacter:
     """
     Calculates the weighted character
 
@@ -660,10 +696,20 @@ def weighted_character(gdf, values, spatial_weights, unique_id, areas=None):
         the name of the left dataframe column, np.array, or pd.Series where is stored area value
 
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    wc : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : GeoDataFrame
+        Series containing used values
+    areas : GeoDataFrame
+        Series containing used areas
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
 
     References
     ----------
@@ -673,61 +719,49 @@ def weighted_character(gdf, values, spatial_weights, unique_id, areas=None):
     Examples
     --------
     >>> sw = libpysal.weights.DistanceBand.from_dataframe(tessellation_df, threshold=100, silence_warnings=True)
-    >>> buildings_df['w_height_100'] = momepy.weighted_character(buildings_df, values='height', spatial_weights=sw,
-                                                                 unique_id='uID')
-    Calculating weighted height...
+    >>> buildings_df['w_height_100'] = momepy.WeightedCharacter(buildings_df, values='height', spatial_weights=sw,
+                                                                 unique_id='uID').wc
     100%|██████████| 144/144 [00:00<00:00, 361.60it/s]
-    Weighted height calculated.
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating weighted character...")
-    data = gdf.copy()
-    if areas is not None:
+    def __init__(self, gdf, values, spatial_weights, unique_id, areas=None):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+
+        data = gdf.copy()
+        if areas is None:
+            areas = gdf.geometry.area
+
         if not isinstance(areas, str):
             data["mm_a"] = areas
             areas = "mm_a"
-    if not isinstance(values, str):
-        data["mm_vals"] = values
-        values = "mm_vals"
-
-    data = data.set_index(unique_id)
+        if not isinstance(values, str):
+            data["mm_vals"] = values
+            values = "mm_vals"
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index]
+        self.areas = data[areas]
+        self.values = data[values]
 
-        if neighbours:
-            building_neighbours = data.loc[neighbours]
+        data = data.set_index(unique_id)
 
-            if areas is not None:
-                results_list.append(
-                    (
-                        sum(building_neighbours[values] * building_neighbours[areas])
-                        + (row[values] * row[areas])
-                    )
-                    / (sum(building_neighbours[areas]) + row[areas])
-                )
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
             else:
-                results_list.append(
-                    (
-                        sum(
-                            building_neighbours[values]
-                            * building_neighbours.geometry.area
-                        )
-                        + (row[values] * row.geometry.area)
-                    )
-                    / (sum(building_neighbours.geometry.area) + row.geometry.area)
-                )
-        else:
-            results_list.append(row[values])
-    series = pd.Series(results_list, index=gdf.index)
+                neighbours = [index]
 
-    print("Weighted character calculated.")
-    return series
+            subset = data.loc[neighbours]
+            results_list.append(
+                (sum(subset[values] * subset[areas])) / (sum(subset[areas]))
+            )
+
+        self.wc = pd.Series(results_list, index=gdf.index)
 
 
-def covered_area(gdf, spatial_weights, unique_id):
+class CoveredArea:
     """
     Calculates the area covered by k steps of morphological tessellation
 
@@ -745,42 +779,47 @@ def covered_area(gdf, spatial_weights, unique_id):
     unique_id : str
         name of the column with unique id used as spatial_weights index.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ca : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
 
     Examples
     --------
     >>> sw = momepy.sw_high(k=3, gdf=tessellation_df, ids='uID')
-    >>> tessellation_df['covered3steps'] = mm.covered_area(tessellation_df, sw, 'uID')
-    Calculating covered area...
+    >>> tessellation_df['covered3steps'] = mm.CoveredArea(tessellation_df, sw, 'uID').ca
     100%|██████████| 144/144 [00:00<00:00, 549.15it/s]
-    Covered area calculated.
 
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating covered area...")
-    data = gdf.set_index(unique_id)
+    def __init__(self, gdf, spatial_weights, unique_id):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = [index]
-        areas = data.loc[neighbours].geometry.area
-        results_list.append(sum(areas))
+        data = gdf.set_index(unique_id)
+
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
 
-    series = pd.Series(results_list, index=gdf.index)
+            areas = data.loc[neighbours].geometry.area
+            results_list.append(sum(areas))
 
-    print("Covered area calculated.")
-    return series
+        self.ca = pd.Series(results_list, index=gdf.index)
 
 
-def wall(gdf, spatial_weights=None):
+class PerimeterWall:
     """
     Calculate the perimeter wall length the joined structure.
 
@@ -792,64 +831,63 @@ def wall(gdf, spatial_weights=None):
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on gdf. It is to denote adjacent buildings (note: based on index, not ID).
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    wall : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw : libpysal.weights
+        spatial weights matrix
 
     Examples
     --------
-    >>> buildings_df['wall_length'] = mm.wall(buildings_df)
-    Calculating perimeter wall length...
+    >>> buildings_df['wall_length'] = mm.PerimeterWall(buildings_df).wall
     Calculating spatial weights...
     Spatial weights ready...
     100%|██████████| 144/144 [00:00<00:00, 4171.39it/s]
-    Perimeter wall length calculated.
 
     Notes
     -----
     It might take a while to compute this character.
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating perimeter wall length...")
+    def __init__(self, gdf, spatial_weights=None):
+        self.gdf = gdf
 
-    # if weights matrix is not passed, generate it from objects
-    if spatial_weights is None:
-        print("Calculating spatial weights...")
-        from libpysal.weights import Queen
+        if spatial_weights is None:
+            print("Calculating spatial weights...")
+            from libpysal.weights import Queen
 
-        spatial_weights = Queen.from_dataframe(gdf, silence_warnings=True)
-        print("Spatial weights ready...")
+            spatial_weights = Queen.from_dataframe(gdf, silence_warnings=True)
+            print("Spatial weights ready...")
+        self.sw = spatial_weights
 
-    # dict to store walls for each uID
-    walls = {}
-    components = pd.Series(spatial_weights.component_labels, index=gdf.index)
+        # dict to store walls for each uID
+        walls = {}
+        components = pd.Series(spatial_weights.component_labels, index=gdf.index)
 
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        # if the id is already present in walls, continue (avoid repetition)
-        if index in walls:
-            continue
-        else:
-            comp = spatial_weights.component_labels[index]
-            to_join = components[components == comp].index
-            joined = gdf.iloc[to_join]
-            dissolved = joined.geometry.buffer(
-                0.01
-            ).unary_union  # buffer to avoid multipolygons where buildings touch by corners only
-            for b in to_join:
-                walls[b] = dissolved.exterior.length  # fill dict with values
-    # copy values from dict to gdf
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        results_list.append(walls[index])
-
-    series = pd.Series(results_list, index=gdf.index)
-    print("Perimeter wall length calculated.")
-    return series
-
-
-def segments_length(gdf, spatial_weights=None, mean=False):
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            # if the id is already present in walls, continue (avoid repetition)
+            if index in walls:
+                continue
+            else:
+                comp = spatial_weights.component_labels[index]
+                to_join = components[components == comp].index
+                joined = gdf.iloc[to_join]
+                dissolved = joined.geometry.buffer(
+                    0.01
+                ).unary_union  # buffer to avoid multipolygons where buildings touch by corners only
+                for b in to_join:
+                    walls[b] = dissolved.exterior.length
+
+        results_list = []
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            results_list.append(walls[index])
+        self.wall = pd.Series(results_list, index=gdf.index)
+
+
+class SegmentsLength:
     """
     Calculate the cummulative or mean length of segments.
 
@@ -867,39 +905,48 @@ def segments_length(gdf, spatial_weights=None, mean=False):
     mean : boolean, optional
         If mean=False it will return total length, if mean=True it will return mean value
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    sl : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw : libpysal.weights
+        spatial weights matrix
+    mean : boolean
+        used mean boolean value
 
     Examples
     --------
-    >>> streets_df['length_neighbours'] = mm.segments_length(streets_df, mean=True)
-    Calculating segments length...
+    >>> streets_df['length_neighbours'] = mm.SegmentsLength(streets_df, mean=True).sl
     Calculating spatial weights...
     Spatial weights ready...
-    Segments length calculated.
     """
-    results_list = []
 
-    print("Calculating segments length...")
+    def __init__(self, gdf, spatial_weights=None, mean=False):
+        self.gdf = gdf
+        self.mean = mean
 
-    if spatial_weights is None:
-        print("Calculating spatial weights...")
-        from libpysal.weights import Queen
+        if spatial_weights is None:
+            print("Calculating spatial weights...")
+            from libpysal.weights import Queen
 
-        spatial_weights = Queen.from_dataframe(gdf)
-        print("Spatial weights ready...")
+            spatial_weights = Queen.from_dataframe(gdf, silence_warnings=True)
+            print("Spatial weights ready...")
+        self.sw = spatial_weights
 
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        neighbours.append(index)
-        dims = gdf.iloc[neighbours].geometry.length
-        if mean:
-            results_list.append(np.mean(dims))
-        else:
-            results_list.append(sum(dims))
+        results_list = []
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
+
+            dims = gdf.iloc[neighbours].geometry.length
+            if mean:
+                results_list.append(np.mean(dims))
+            else:
+                results_list.append(sum(dims))
 
-    series = pd.Series(results_list, index=gdf.index)
-    print("Segments length calculated.")
-    return series
+        self.sl = pd.Series(results_list, index=gdf.index)
diff --git a/momepy/distribution.py b/momepy/distribution.py
index 27ee44f..3d241e0 100644
--- a/momepy/distribution.py
+++ b/momepy/distribution.py
@@ -11,7 +11,7 @@ import pandas as pd
 import statistics
 
 
-def orientation(gdf):
+class Orientation:
     """
     Calculate the orientation of object
 
@@ -23,10 +23,12 @@ def orientation(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects to analyse
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    o : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     References
     ---------
@@ -35,78 +37,75 @@ def orientation(gdf):
 
     Examples
     --------
-    >>> buildings_df['orientation'] = momepy.orientation(buildings_df)
-    Calculating orientations...
+    >>> buildings_df['orientation'] = momepy.Orientation(buildings_df).o
     100%|██████████| 144/144 [00:00<00:00, 630.54it/s]
-    Orientations calculated.
     >>> buildings_df['orientation'][0]
     41.05146788287027
     """
-    # define empty list for results
-    results_list = []
-
-    print("Calculating orientations...")
-
-    def _azimuth(point1, point2):
-        """azimuth between 2 shapely points (interval 0 - 180)"""
-        angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
-        return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
-
-    # iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        bbox = list(row["geometry"].minimum_rotated_rectangle.exterior.coords)
-        centroid_ab = LineString([bbox[0], bbox[1]]).centroid
-        centroid_cd = LineString([bbox[2], bbox[3]]).centroid
-        axis1 = centroid_ab.distance(centroid_cd)
-
-        centroid_bc = LineString([bbox[1], bbox[2]]).centroid
-        centroid_da = LineString([bbox[3], bbox[0]]).centroid
-        axis2 = centroid_bc.distance(centroid_da)
-
-        if axis1 <= axis2:
-            az = _azimuth(centroid_bc, centroid_da)
-            if 90 > az >= 45:
-                diff = az - 45
-                az = az - 2 * diff
-            elif 135 > az >= 90:
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            elif 181 > az >= 135:
-                diff = az - 135
-                az = az - 2 * diff
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            results_list.append(az)
-        else:
-            az = 170
-            az = _azimuth(centroid_ab, centroid_cd)
-            if 90 > az >= 45:
-                diff = az - 45
-                az = az - 2 * diff
-            elif 135 > az >= 90:
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            elif 181 > az >= 135:
-                diff = az - 135
-                az = az - 2 * diff
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            results_list.append(az)
-
-    series = pd.Series(results_list, index=gdf.index)
-    print("Orientations calculated.")
-    return series
 
-
-def shared_walls_ratio(gdf, unique_id, perimeters=None):
+    def __init__(self, gdf):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
+
+        def _azimuth(point1, point2):
+            """azimuth between 2 shapely points (interval 0 - 180)"""
+            angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
+            return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
+
+        # iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            bbox = list(row["geometry"].minimum_rotated_rectangle.exterior.coords)
+            centroid_ab = LineString([bbox[0], bbox[1]]).centroid
+            centroid_cd = LineString([bbox[2], bbox[3]]).centroid
+            axis1 = centroid_ab.distance(centroid_cd)
+
+            centroid_bc = LineString([bbox[1], bbox[2]]).centroid
+            centroid_da = LineString([bbox[3], bbox[0]]).centroid
+            axis2 = centroid_bc.distance(centroid_da)
+
+            if axis1 <= axis2:
+                az = _azimuth(centroid_bc, centroid_da)
+                if 90 > az >= 45:
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 135 > az >= 90:
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 181 > az >= 135:
+                    diff = az - 135
+                    az = az - 2 * diff
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                results_list.append(az)
+            else:
+                az = 170
+                az = _azimuth(centroid_ab, centroid_cd)
+                if 90 > az >= 45:
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 135 > az >= 90:
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 181 > az >= 135:
+                    diff = az - 135
+                    az = az - 2 * diff
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                results_list.append(az)
+
+        self.o = pd.Series(results_list, index=gdf.index)
+
+
+class SharedWallsRatio:
     """
     Calculate shared walls ratio
 
@@ -122,10 +121,18 @@ def shared_walls_ratio(gdf, unique_id, perimeters=None):
     perimeters : str, list, np.array, pd.Series (default None)
         the name of the dataframe column, np.array, or pd.Series where is stored perimeter value
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    swr : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    id : Series
+        Series containing used unique ID
+    perimeters : GeoDataFrame
+        Series containing used perimeters values
+    sindex : rtree spatial index
+        spatial index of gdf
 
     References
     ---------
@@ -136,60 +143,53 @@ def shared_walls_ratio(gdf, unique_id, perimeters=None):
 
     Examples
     --------
-    >>> buildings_df['swr'] = momepy.shared_walls_ratio(buildings_df, 'uID')
-    Generating spatial index...
-    Calculating shared walls ratio...
+    >>> buildings_df['swr'] = momepy.SharedWallsRatio(buildings_df, 'uID').swr
     100%|██████████| 144/144 [00:00<00:00, 648.72it/s]
-    Shared walls ratio calculated.
     >>> buildings_df['swr'][10]
     0.3424804411228673
     """
-    print("Generating spatial index...")
-    gdf = gdf.copy()
-    sindex = gdf.sindex  # define rtree index
-    # define empty list for results
-    results_list = []
-
-    print("Calculating shared walls ratio...")
-
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+
+    def __init__(self, gdf, unique_id, perimeters=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        self.sindex = gdf.sindex  # define rtree index
+        # define empty list for results
+        results_list = []
+
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
-    if not isinstance(unique_id, str):
-        gdf["mm_uid"] = unique_id
-        unique_id = "mm_uid"
-
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        neighbors = list(sindex.intersection(row.geometry.bounds))
-        neighbors.remove(index)
-
-        # if no neighbour exists
-        length = 0
-        if not neighbors:
-            results_list.append(0)
         else:
-            for i in neighbors:
-                subset = gdf.iloc[i]["geometry"]
-                length = length + row.geometry.intersection(subset).length
-            results_list.append(length / row[perimeters])
-    series = pd.Series(results_list, index=gdf.index)
-    print("Shared walls ratio calculated.")
-
-    return series
-
-
-def street_alignment(
-    left,
-    right,
-    orientations,
-    network_id=None,
-    left_network_id=None,
-    right_network_id=None,
-):
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+
+        self.perimeters = gdf[perimeters]
+
+        if not isinstance(unique_id, str):
+            gdf["mm_uid"] = unique_id
+            unique_id = "mm_uid"
+        self.id = gdf[unique_id]
+
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            neighbors = list(self.sindex.intersection(row.geometry.bounds))
+            neighbors.remove(index)
+
+            # if no neighbour exists
+            length = 0
+            if not neighbors:
+                results_list.append(0)
+            else:
+                for i in neighbors:
+                    subset = gdf.iloc[i]["geometry"]
+                    length = length + row.geometry.intersection(subset).length
+                results_list.append(length / row[perimeters])
+
+        self.swr = pd.Series(results_list, index=gdf.index)
+
+
+class StreetAlignment:
     """
     Calculate the difference between street orientation and orientation of object in degrees
 
@@ -218,90 +218,110 @@ def street_alignment(
         the name of the right dataframe column, np.array, or pd.Series of streets with unique network id (has to be defined beforehand)
         (can be defined using :py:func:`momepy.elements.unique_id`)
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    sa : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    network_id : str
+        the name of the column storing network ID in both left and right
+    left_network_id : Series
+        Series containing used left ID
+    right_network_id : Series
+        Series containing used right ID
 
     Examples
     --------
-    >>> buildings_df['street_alignment'] = momepy.street_alignment(buildings_df, streets_df, 'orientation', 'nID', 'nID')
-    Calculating street alignments...
+    >>> buildings_df['street_alignment'] = momepy.StreetAlignment(buildings_df, streets_df, 'orientation', 'nID', 'nID').sa
     100%|██████████| 144/144 [00:00<00:00, 529.94it/s]
-    Street alignments calculated.
     >>> buildings_df['street_alignment'][0]
     0.29073888476702336
     """
-    # define empty list for results
-    results_list = []
-
-    print("Calculating street alignments...")
-    left = left.copy()
-    right = right.copy()
-
-    if network_id:
-        left_network_id = network_id
-        right_network_id = network_id
-    else:
-        if left_network_id is None and right_network_id is not None:
-            raise ValueError("left_network_id not set.")
-        if left_network_id is not None and right_network_id is None:
-            raise ValueError("right_network_id not set.")
-        if left_network_id is None and right_network_id is None:
-            raise ValueError(
-                "Network ID not set. Use either network_id or left_network_id and right_network_id."
-            )
 
-    if not isinstance(orientations, str):
-        left["mm_o"] = orientations
-        orientations = "mm_o"
-    if not isinstance(left_network_id, str):
-        left["mm_nid"] = left_network_id
-        left_network_id = "mm_nid"
-    if not isinstance(right_network_id, str):
-        right["mm_nis"] = right_network_id
-        right_network_id = "mm_nis"
-
-    def azimuth(point1, point2):
-        """azimuth between 2 shapely points (interval 0 - 180)"""
-        angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
-        return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
-
-    # iterating over rows one by one
-    for index, row in tqdm(left.iterrows(), total=left.shape[0]):
-        if pd.isnull(row[left_network_id]):
-            results_list.append(0)
+    def __init__(
+        self,
+        left,
+        right,
+        orientations,
+        network_id=None,
+        left_network_id=None,
+        right_network_id=None,
+    ):
+        self.left = left
+        self.right = right
+        self.network_id = network_id
+
+        # define empty list for results
+        results_list = []
+
+        left = left.copy()
+        right = right.copy()
+
+        if network_id:
+            left_network_id = network_id
+            right_network_id = network_id
         else:
-            network_id = row[left_network_id]
-            streetssub = right.loc[right[right_network_id] == network_id]
-            start = Point(streetssub.iloc[0]["geometry"].coords[0])
-            end = Point(streetssub.iloc[0]["geometry"].coords[-1])
-            az = azimuth(start, end)
-            if 90 > az >= 45:
-                diff = az - 45
-                az = az - 2 * diff
-            elif 135 > az >= 90:
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            elif 181 > az >= 135:
-                diff = az - 135
-                az = az - 2 * diff
-                diff = az - 90
-                az = az - 2 * diff
-                diff = az - 45
-                az = az - 2 * diff
-            results_list.append(abs(row[orientations] - az))
-    series = pd.Series(results_list, index=left.index)
-
-    print("Street alignments calculated.")
-    return series
-
-
-def cell_alignment(
-    left, right, left_orientations, right_orientations, left_unique_id, right_unique_id
-):
+            if left_network_id is None and right_network_id is not None:
+                raise ValueError("left_network_id not set.")
+            if left_network_id is not None and right_network_id is None:
+                raise ValueError("right_network_id not set.")
+            if left_network_id is None and right_network_id is None:
+                raise ValueError(
+                    "Network ID not set. Use either network_id or left_network_id and right_network_id."
+                )
+
+        if not isinstance(orientations, str):
+            left["mm_o"] = orientations
+            orientations = "mm_o"
+        self.orientations = left[orientations]
+
+        if not isinstance(left_network_id, str):
+            left["mm_nid"] = left_network_id
+            left_network_id = "mm_nid"
+        self.left_network_id = left[left_network_id]
+        if not isinstance(right_network_id, str):
+            right["mm_nis"] = right_network_id
+            right_network_id = "mm_nis"
+        self.right_network_id = right[right_network_id]
+
+        def _azimuth(point1, point2):
+            """azimuth between 2 shapely points (interval 0 - 180)"""
+            angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
+            return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
+
+        # iterating over rows one by one
+        for index, row in tqdm(left.iterrows(), total=left.shape[0]):
+            if pd.isnull(row[left_network_id]):
+                results_list.append(0)
+            else:
+                network_id = row[left_network_id]
+                streetssub = right.loc[right[right_network_id] == network_id]
+                start = Point(streetssub.iloc[0]["geometry"].coords[0])
+                end = Point(streetssub.iloc[0]["geometry"].coords[-1])
+                az = _azimuth(start, end)
+                if 90 > az >= 45:
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 135 > az >= 90:
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                elif 181 > az >= 135:
+                    diff = az - 135
+                    az = az - 2 * diff
+                    diff = az - 90
+                    az = az - 2 * diff
+                    diff = az - 45
+                    az = az - 2 * diff
+                results_list.append(abs(row[orientations] - az))
+        self.sa = pd.Series(results_list, index=left.index)
+
+
+class CellAlignment:
     """
     Calculate the difference between cell orientation and orientation of object
 
@@ -325,52 +345,75 @@ def cell_alignment(
     right_unique_id : str
         the name of the right dataframe column with unique id shared between left and right gdf
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ca : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    left_orientations : Series
+        Series containing used left orientations
+    right_orientations : Series
+        Series containing used right orientations
+    left_unique_id : Series
+        Series containing used left ID
+    right_unique_id : Series
+        Series containing used right ID
 
     Examples
     --------
-    >>> buildings_df['cell_alignment'] = momepy.cell_alignment(buildings_df, tessellation_df, 'bl_orient', 'tes_orient', 'uID', 'uID')
-    Calculating cell alignments...
+    >>> buildings_df['cell_alignment'] = momepy.CellAlignment(buildings_df, tessellation_df, 'bl_orient', 'tes_orient', 'uID', 'uID').ca
     100%|██████████| 144/144 [00:00<00:00, 799.09it/s]
-    Cell alignments calculated.
     >>> buildings_df['cell_alignment'][0]
     0.8795123936951939
 
     """
-    print("Calculating cell alignments...")
-    left = left.copy()
-    right = right.copy()
-    if not isinstance(left_orientations, str):
-        left["mm_o"] = left_orientations
-        left_orientations = "mm_o"
-    if not isinstance(right_orientations, str):
-        right["mm_o"] = right_orientations
-        right_orientations = "mm_o"
-
-    # define empty list for results
-    results_list = []
-
-    for index, row in tqdm(left.iterrows(), total=left.shape[0]):
-
-        results_list.append(
-            abs(
-                row[left_orientations]
-                - right[right[right_unique_id] == row[left_unique_id]][
-                    right_orientations
-                ].iloc[0]
-            )
-        )
 
-    series = pd.Series(results_list, index=left.index)
+    def __init__(
+        self,
+        left,
+        right,
+        left_orientations,
+        right_orientations,
+        left_unique_id,
+        right_unique_id,
+    ):
+        self.left = left
+        self.right = right
+
+        left = left.copy()
+        right = right.copy()
+        if not isinstance(left_orientations, str):
+            left["mm_o"] = left_orientations
+            left_orientations = "mm_o"
+        self.left_orientations = left[left_orientations]
+        if not isinstance(right_orientations, str):
+            right["mm_o"] = right_orientations
+            right_orientations = "mm_o"
+        self.right_orientations = right[right_orientations]
+        self.left_unique_id = left[left_unique_id]
+        self.right_unique_id = right[right_unique_id]
+
+        # define empty list for results
+        results_list = []
+
+        for index, row in tqdm(left.iterrows(), total=left.shape[0]):
 
-    print("Cell alignments calculated.")
-    return series
+            results_list.append(
+                abs(
+                    row[left_orientations]
+                    - right[right[right_unique_id] == row[left_unique_id]][
+                        right_orientations
+                    ].iloc[0]
+                )
+            )
 
+        self.ca = pd.Series(results_list, index=left.index)
 
-def alignment(gdf, spatial_weights, unique_id, orientations):
+
+class Alignment:
 
     """
     Calculate the mean deviation of solar orientation of objects on adjacent cells from an object
@@ -390,51 +433,58 @@ def alignment(gdf, spatial_weights, unique_id, orientations):
     unique_id : str
         name of the column with unique id used as spatial_weights index.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    a : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    orientations : Series
+        Series containing used orientation values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
 
     Examples
     --------
-    >>> buildings_df['alignment'] = momepy.alignment(buildings_df, sw, 'uID', bl_orient)
-    Calculating alignments...
-    Calculating spatial weights...
-    Spatial weights ready...
+    >>> buildings_df['alignment'] = momepy.Alignment(buildings_df, sw, 'uID', bl_orient).a
     100%|██████████| 144/144 [00:01<00:00, 140.84it/s]
-    Alignments calculated.
     >>> buildings_df['alignment'][0]
     18.299481296455237
     """
-    # define empty list for results
-    results_list = []
-    gdf = gdf.copy()
-    if not isinstance(orientations, str):
-        gdf["mm_o"] = orientations
-        orientations = "mm_o"
 
-    print("Calculating alignments...")
-    data = gdf.set_index(unique_id)
+    def __init__(self, gdf, spatial_weights, unique_id, orientations):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
 
-    # iterating over rows one by one
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+        # define empty list for results
+        results_list = []
+        gdf = gdf.copy()
+        if not isinstance(orientations, str):
+            gdf["mm_o"] = orientations
+            orientations = "mm_o"
+        self.orientations = gdf[orientations]
 
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            orientation = data.loc[neighbours][orientations]
-            deviations = abs(orientation - row[orientations])
+        data = gdf.set_index(unique_id)
 
-            results_list.append(statistics.mean(deviations))
-        else:
-            results_list.append(0)
+        # iterating over rows one by one
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
 
-    series = pd.Series(results_list, index=gdf.index)
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                orientation = data.loc[neighbours][orientations]
+                deviations = abs(orientation - row[orientations])
 
-    print("Alignments calculated.")
-    return series
+                results_list.append(statistics.mean(deviations))
+            else:
+                results_list.append(0)
 
+        self.a = pd.Series(results_list, index=gdf.index)
 
-def neighbour_distance(gdf, spatial_weights, unique_id):
+
+class NeighborDistance:
     """
     Calculate the mean distance to adjacent buildings (based on spatial_weights)
 
@@ -450,10 +500,16 @@ def neighbour_distance(gdf, spatial_weights, unique_id):
     unique_id : str
         name of the column with unique id used as spatial_weights index.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    nd : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
 
     References
     ---------
@@ -462,39 +518,36 @@ def neighbour_distance(gdf, spatial_weights, unique_id):
 
     Examples
     --------
-    >>> buildings_df['neighbour_distance'] = momepy.neighbour_distance(buildings_df, sw, 'uID')
-    Calculating distances...
+    >>> buildings_df['neighbour_distance'] = momepy.NeighborDistance(buildings_df, sw, 'uID').nd
     100%|██████████| 144/144 [00:00<00:00, 345.78it/s]
-    Distances calculated.
     >>> buildings_df['neighbour_distance'][0]
     29.18589019096464
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating distances...")
-    data = gdf.set_index(unique_id)
+    def __init__(self, gdf, spatial_weights, unique_id):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        # define empty list for results
+        results_list = []
 
-    # iterating over rows one by one
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index]
-        building_neighbours = data.loc[neighbours]
-        if len(building_neighbours) > 0:
-            results_list.append(
-                np.mean(building_neighbours.geometry.distance(row["geometry"]))
-            )
-        else:
-            results_list.append(0)
+        data = gdf.set_index(unique_id)
 
-    series = pd.Series(results_list, index=gdf.index)
+        # iterating over rows one by one
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index]
+            building_neighbours = data.loc[neighbours]
+            if len(building_neighbours) > 0:
+                results_list.append(
+                    np.mean(building_neighbours.geometry.distance(row["geometry"]))
+                )
+            else:
+                results_list.append(0)
 
-    print("Distances calculated.")
-    return series
+        self.nd = pd.Series(results_list, index=gdf.index)
 
 
-def mean_interbuilding_distance(
-    gdf, spatial_weights, unique_id, spatial_weights_higher=None, order=3
-):
+class MeanInterbuildingDistance:
     """
     Calculate the mean interbuilding distance within x topological steps
 
@@ -518,10 +571,20 @@ def mean_interbuilding_distance(
     order : int
         Order of Queen contiguity
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    mid : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    sw_higher : libpysal.weights
+        Spatial weights matrix of higher order
+    order : int
+        Order of Queen contiguity (only if spatial_weights_higher was not set)
 
     Notes
     -----
@@ -529,73 +592,79 @@ def mean_interbuilding_distance(
 
     Examples
     --------
-    >>> buildings_df['mean_interbuilding_distance'] = momepy.mean_interbuilding_distance(buildings_df, sw, 'uID')
-    Calculating mean interbuilding distances...
+    >>> buildings_df['mean_interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings_df, sw, 'uID').mid
     Generating weights matrix (Queen) of 3 topological steps...
     Generating adjacency matrix based on weights matrix...
     Computing interbuilding distances...
     100%|██████████| 746/746 [00:03<00:00, 200.14it/s]
     Computing mean interbuilding distances...
     100%|██████████| 144/144 [00:00<00:00, 317.42it/s]
-    Mean interbuilding distances calculated.
     >>> buildings_df['mean_interbuilding_distance'][0]
     29.305457092042744
     """
 
-    print("Calculating mean interbuilding distances...")
-    if spatial_weights_higher is None:
-        print(
-            "Generating weights matrix (Queen) of {} topological steps...".format(order)
-        )
-        from momepy import sw_high
-
-        # matrix to define area of analysis (more steps)
-        spatial_weights_higher = sw_high(k=order, weights=spatial_weights)
-
-    # define empty list for results
-    results_list = []
-
-    print("Generating adjacency matrix based on weights matrix...")
-    # define adjacency list from lipysal
-    adj_list = spatial_weights.to_adjlist()
-    adj_list["distance"] = -1
-
-    print("Computing interbuilding distances...")
-    # measure each interbuilding distance of neighbours and save them to adjacency list
-    for index, row in tqdm(adj_list.iterrows(), total=adj_list.shape[0]):
-        inverted = adj_list[(adj_list.focal == row.neighbor)][
-            (adj_list.neighbor == row.focal)
-        ].iloc[0]["distance"]
-        if inverted == -1:
-            building_object = gdf.loc[gdf[unique_id] == row.focal.astype(int)]
-
-            building_neighbour = gdf.loc[gdf[unique_id] == row.neighbor.astype(int)]
-            adj_list.loc[index, "distance"] = building_neighbour.iloc[
-                0
-            ].geometry.distance(building_object.iloc[0].geometry)
-        else:
-            adj_list.at[index, "distance"] = inverted
-
-    print("Computing mean interbuilding distances...")
-    # iterate over objects to get the final values
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        # id to match spatial weights
-        uid = row[unique_id]
-        # define neighbours based on weights matrix defining analysis area
-        neighbours = spatial_weights_higher.neighbors[uid].copy()
-        neighbours.append(uid)
-        if neighbours:
-            selection = adj_list[adj_list.focal.isin(neighbours)][
-                adj_list.neighbor.isin(neighbours)
-            ]
-            results_list.append(np.nanmean(selection.distance))
-
-    series = pd.Series(results_list, index=gdf.index)
-    print("Mean interbuilding distances calculated.")
-    return series
-
-
-def neighbouring_street_orientation_deviation(gdf):
+    def __init__(
+        self, gdf, spatial_weights, unique_id, spatial_weights_higher=None, order=3
+    ):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+
+        if spatial_weights_higher is None:
+            print(
+                "Generating weights matrix (Queen) of {} topological steps...".format(
+                    order
+                )
+            )
+            self.order = order
+            from momepy import sw_high
+
+            # matrix to define area of analysis (more steps)
+            spatial_weights_higher = sw_high(k=order, weights=spatial_weights)
+        self.sw_higher = spatial_weights_higher
+
+        # define empty list for results
+        results_list = []
+
+        print("Generating adjacency matrix based on weights matrix...")
+        # define adjacency list from lipysal
+        adj_list = spatial_weights.to_adjlist()
+        adj_list["distance"] = -1
+
+        print("Computing interbuilding distances...")
+        # measure each interbuilding distance of neighbours and save them to adjacency list
+        for index, row in tqdm(adj_list.iterrows(), total=adj_list.shape[0]):
+            inverted = adj_list[(adj_list.focal == row.neighbor)][
+                (adj_list.neighbor == row.focal)
+            ].iloc[0]["distance"]
+            if inverted == -1:
+                building_object = gdf.loc[gdf[unique_id] == row.focal.astype(int)]
+
+                building_neighbour = gdf.loc[gdf[unique_id] == row.neighbor.astype(int)]
+                adj_list.loc[index, "distance"] = building_neighbour.iloc[
+                    0
+                ].geometry.distance(building_object.iloc[0].geometry)
+            else:
+                adj_list.at[index, "distance"] = inverted
+
+        print("Computing mean interbuilding distances...")
+        # iterate over objects to get the final values
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            # id to match spatial weights
+            uid = row[unique_id]
+            # define neighbours based on weights matrix defining analysis area
+            neighbours = spatial_weights_higher.neighbors[uid].copy()
+            neighbours.append(uid)
+            if neighbours:
+                selection = adj_list[adj_list.focal.isin(neighbours)][
+                    adj_list.neighbor.isin(neighbours)
+                ]
+                results_list.append(np.nanmean(selection.distance))
+
+        self.mid = pd.Series(results_list, index=gdf.index)
+
+
+class NeighboringStreetOrientationDeviation:
     """
     Calculate the mean deviation of solar orientation of adjacent streets
 
@@ -610,89 +679,95 @@ def neighbouring_street_orientation_deviation(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing street network to analyse
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    nsod : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    orientation : Series
+        Series containing used street orientation values
+    sindex : rtree spatial index
+        spatial index of gdf
 
     Examples
     --------
-    >>> streets_df['orient_dev'] = momepy.neighbouring_street_orientation_deviation(streets_df)
-    Calculating street alignments...
+    >>> streets_df['orient_dev'] = momepy.NeighboringStreetOrientationDeviation(streets_df).nsod
     Preparing street orientations...
     Generating spatial index...
     100%|██████████| 33/33 [00:00<00:00, 249.02it/s]
-    Street alignments calculated.
     >>> streets_df['orient_dev'][6]
     7.043096518688273
     """
-    # define empty list for results
-    results_list = []
-    gdf = gdf.copy()
-    print("Calculating street alignments...")
-
-    def azimuth(point1, point2):
-        """azimuth between 2 shapely points (interval 0 - 180)"""
-        angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
-        return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
-
-    # iterating over rows one by one
-    print(" Preparing street orientations...")
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-
-        start = Point(row["geometry"].coords[0])
-        end = Point(row["geometry"].coords[-1])
-        az = azimuth(start, end)
-        if 90 > az >= 45:
-            diff = az - 45
-            az = az - 2 * diff
-        elif 135 > az >= 90:
-            diff = az - 90
-            az = az - 2 * diff
-            diff = az - 45
-            az = az - 2 * diff
-        elif 181 > az >= 135:
-            diff = az - 135
-            az = az - 2 * diff
-            diff = az - 90
-            az = az - 2 * diff
-            diff = az - 45
-            az = az - 2 * diff
-        results_list.append(az)
-    series = pd.Series(results_list, index=gdf.index)
-
-    gdf["tmporient"] = series
-
-    print(" Generating spatial index...")
-    sindex = gdf.sindex
-    results_list = []
-
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        possible_neighbors_idx = list(sindex.intersection(row.geometry.bounds))
-        possible_neighbours = gdf.iloc[possible_neighbors_idx]
-        neighbors = possible_neighbours[possible_neighbours.intersects(row.geometry)]
-        neighbors.drop([index])
-
-        orientations = []
-        for idx, r in neighbors.iterrows():
-            orientations.append(r.tmporient)
-
-        deviations = []
-        for o in orientations:
-            dev = abs(o - row.tmporient)
-            deviations.append(dev)
-
-        if deviations:
-            results_list.append(np.mean(deviations))
-        else:
-            results_list.append(0)
 
-    series = pd.Series(results_list, index=gdf.index)
-    print("Street alignments calculated.")
-    return series
+    def __init__(self, gdf):
+        self.gdf = gdf
+
+        results_list = []
+        gdf = gdf.copy()
+
+        def _azimuth(point1, point2):
+            """azimuth between 2 shapely points (interval 0 - 180)"""
+            angle = np.arctan2(point2.x - point1.x, point2.y - point1.y)
+            return np.degrees(angle) if angle > 0 else np.degrees(angle) + 180
+
+        # iterating over rows one by one
+        print(" Preparing street orientations...")
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
 
+            start = Point(row["geometry"].coords[0])
+            end = Point(row["geometry"].coords[-1])
+            az = _azimuth(start, end)
+            if 90 > az >= 45:
+                diff = az - 45
+                az = az - 2 * diff
+            elif 135 > az >= 90:
+                diff = az - 90
+                az = az - 2 * diff
+                diff = az - 45
+                az = az - 2 * diff
+            elif 181 > az >= 135:
+                diff = az - 135
+                az = az - 2 * diff
+                diff = az - 90
+                az = az - 2 * diff
+                diff = az - 45
+                az = az - 2 * diff
+            results_list.append(az)
+        self.orientation = pd.Series(results_list, index=gdf.index)
+
+        gdf["tmporient"] = self.orientation
+
+        print(" Generating spatial index...")
+        self.sindex = gdf.sindex
+        results_list = []
+
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            possible_neighbors_idx = list(self.sindex.intersection(row.geometry.bounds))
+            possible_neighbours = gdf.iloc[possible_neighbors_idx]
+            neighbors = possible_neighbours[
+                possible_neighbours.intersects(row.geometry)
+            ]
+            neighbors.drop([index])
+
+            orientations = []
+            for idx, r in neighbors.iterrows():
+                orientations.append(r.tmporient)
+
+            deviations = []
+            for o in orientations:
+                dev = abs(o - row.tmporient)
+                deviations.append(dev)
+
+            if deviations:
+                results_list.append(np.mean(deviations))
+            else:
+                results_list.append(0)
+
+        self.nsod = pd.Series(results_list, index=gdf.index)
 
-def building_adjacency(gdf, spatial_weights_higher, unique_id, spatial_weights=None):
+
+class BuildingAdjacency:
     """
     Calculate the level of building adjacency
 
@@ -715,10 +790,18 @@ def building_adjacency(gdf, spatial_weights_higher, unique_id, spatial_weights=N
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on gdf. It is to denote adjacent buildings (note: based on unique ID).
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ba : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    sw_higher : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    sw : libpysal.weights
+        spatial weights matrix
 
     References
     ---------
@@ -727,55 +810,53 @@ def building_adjacency(gdf, spatial_weights_higher, unique_id, spatial_weights=N
 
     Examples
     --------
-    >>> buildings_df['adjacency'] = momepy.building_adjacency(buildings_df, swh, unique_id='uID')
-    Calculating adjacency...
+    >>> buildings_df['adjacency'] = momepy.BuildingAdjacency(buildings_df, swh, unique_id='uID').ba
     Calculating spatial weights...
     Spatial weights ready...
     Generating weights matrix (Queen) of 3 topological steps...
     100%|██████████| 144/144 [00:00<00:00, 9301.73it/s]
     Calculating adjacency within k steps...
     100%|██████████| 144/144 [00:00<00:00, 335.55it/s]
-    Adjacency calculated.
     >>> buildings_df['adjacency'][10]
     0.23809523809523808
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating adjacency...")
+    def __init__(self, gdf, spatial_weights_higher, unique_id, spatial_weights=None):
+        self.gdf = gdf
+        self.sw_higher = spatial_weights_higher
+        self.id = gdf[unique_id]
+        results_list = []
 
-    # if weights matrix is not passed, generate it from gdf
-    if spatial_weights is None:
-        print("Calculating spatial weights...")
-        from libpysal.weights import Queen
+        # if weights matrix is not passed, generate it from gdf
+        if spatial_weights is None:
+            print("Calculating spatial weights...")
+            from libpysal.weights import Queen
 
-        spatial_weights = Queen.from_dataframe(
-            gdf, silence_warnings=True, ids=unique_id
-        )
-        print("Spatial weights ready...")
-
-    patches = dict(zip(gdf[unique_id], spatial_weights.component_labels))
+            spatial_weights = Queen.from_dataframe(
+                gdf, silence_warnings=True, ids=unique_id
+            )
+            print("Spatial weights ready...")
 
-    print("Calculating adjacency within k steps...")
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        neighbours = spatial_weights_higher.neighbors[row[unique_id]].copy()
-        if neighbours:
-            neighbours.append(row[unique_id])
+        self.sw = spatial_weights
+        patches = dict(zip(gdf[unique_id], spatial_weights.component_labels))
 
-            patches_sub = [patches[x] for x in neighbours]
-            patches_nr = len(set(patches_sub))
+        print("Calculating adjacency within k steps...")
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            neighbours = spatial_weights_higher.neighbors[row[unique_id]].copy()
+            if neighbours:
+                neighbours.append(row[unique_id])
 
-            results_list.append(patches_nr / len(neighbours))
-        else:
-            results_list.append(0)
+                patches_sub = [patches[x] for x in neighbours]
+                patches_nr = len(set(patches_sub))
 
-    series = pd.Series(results_list, index=gdf.index)
+                results_list.append(patches_nr / len(neighbours))
+            else:
+                results_list.append(0)
 
-    print("Adjacency calculated.")
-    return series
+        self.ba = pd.Series(results_list, index=gdf.index)
 
 
-def neighbours(gdf, spatial_weights, unique_id, weighted=False):
+class Neighbors:
     """
     Calculate the number of topological neighbours of each object.
 
@@ -796,10 +877,20 @@ def neighbours(gdf, spatial_weights, unique_id, weighted=False):
     weighted : bool (default False)
         if weighted=True, number of neighbours will be divided by the perimeter of object, to return relative value
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    n : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    weighted : bool
+        used weighted value
 
     References
     ---------
@@ -810,24 +901,25 @@ def neighbours(gdf, spatial_weights, unique_id, weighted=False):
     Examples
     --------
     >>> sw = libpysal.weights.contiguity.Queen.from_dataframe(tessellation_df, ids='uID')
-    >>> tessellation_df['neighbours'] = momepy.neighbours(tessellation_df, sw, 'uID')
-    Calculating neighbours...
+    >>> tessellation_df['neighbours'] = momepy.Neighbors(tessellation_df, sw, 'uID').n
     100%|██████████| 144/144 [00:00<00:00, 6909.50it/s]
-    Neighbours calculated.
     >>> tessellation_df['neighbours'][0]
     4
     """
 
-    print("Calculating neighbours...")
-    neighbours = []
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        if weighted is True:
-            neighbours.append(
-                spatial_weights.cardinalities[row[unique_id]] / row.geometry.length
-            )
-        else:
-            neighbours.append(spatial_weights.cardinalities[row[unique_id]])
-
-    series = pd.Series(neighbours, index=gdf.index)
-    print("Neighbours calculated.")
-    return series
+    def __init__(self, gdf, spatial_weights, unique_id, weighted=False):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.weighted = weighted
+
+        neighbours = []
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            if weighted is True:
+                neighbours.append(
+                    spatial_weights.cardinalities[row[unique_id]] / row.geometry.length
+                )
+            else:
+                neighbours.append(spatial_weights.cardinalities[row[unique_id]])
+
+        self.n = pd.Series(neighbours, index=gdf.index)
diff --git a/momepy/diversity.py b/momepy/diversity.py
index b83559e..89bc732 100644
--- a/momepy/diversity.py
+++ b/momepy/diversity.py
@@ -10,7 +10,7 @@ import scipy as sp
 import numpy as np
 
 
-def rng(gdf, values, spatial_weights, unique_id, rng=(0, 100), **kwargs):
+class Range:
     """
     Calculates the range of values within neighbours defined in `spatial_weights`.
 
@@ -35,10 +35,22 @@ def rng(gdf, values, spatial_weights, unique_id, rng=(0, 100), **kwargs):
     **kwargs : keyword arguments
         optional arguments for `scipy.stats.iqr`
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    r : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    rng : tuple
+        range
+    kwargs : dict
+        kwargs
 
     References
     ----------
@@ -48,42 +60,43 @@ def rng(gdf, values, spatial_weights, unique_id, rng=(0, 100), **kwargs):
     Examples
     --------
     >>> sw = momepy.sw_high(k=3, gdf=tessellation_df, ids='uID')
-    >>> tessellation_df['area_IQR_3steps'] = mm.rng(tessellation_df, 'area', sw, 'uID', rng=(25, 75))
-    Calculating range...
+    >>> tessellation_df['area_IQR_3steps'] = mm.Range(tessellation_df, 'area', sw, 'uID', rng=(25, 75)).r
     100%|██████████| 144/144 [00:00<00:00, 722.50it/s]
-    Range calculated.
 
 
     """
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-    print("Calculating range...")
 
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
+    def __init__(self, gdf, values, spatial_weights, unique_id, rng=(0, 100), **kwargs):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.rng = rng
+        self.kwargs = kwargs
 
-    data = data.set_index(unique_id)
+        data = gdf.copy()
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = [index]
-        values_list = data.loc[neighbours][values]
+        data = data.set_index(unique_id)
 
-        results_list.append(sp.stats.iqr(values_list, rng=rng, **kwargs))
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
 
-    series = pd.Series(results_list, index=gdf.index)
+            values_list = data.loc[neighbours][values]
+            results_list.append(sp.stats.iqr(values_list, rng=rng, **kwargs))
 
-    print("Range calculated.")
-    return series
+        self.r = pd.Series(results_list, index=gdf.index)
 
 
-def theil(gdf, values, spatial_weights, unique_id, rng=None):
+class Theil:
     """
     Calculates the Theil measure of inequality of values within neighbours defined in `spatial_weights`.
 
@@ -106,71 +119,71 @@ def theil(gdf, values, spatial_weights, unique_id, rng=None):
         Percentiles over which to compute the range. Each must be
         between 0 and 100, inclusive. The order of the elements is not important.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    t : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    rng : tuple
+        range
 
     Examples
     --------
     >>> sw = momepy.sw_high(k=3, gdf=tessellation_df, ids='uID')
-    >>> tessellation_df['area_Theil'] = mm.theil(tessellation_df, 'area', sw, 'uID')
-    Calculating Theil index...
+    >>> tessellation_df['area_Theil'] = mm.Theil(tessellation_df, 'area', sw, 'uID').t
     100%|██████████| 144/144 [00:00<00:00, 597.37it/s]
-    Theil index calculated.
     """
-    try:
-        from inequality.theil import Theil
-    except ImportError:
+
+    def __init__(self, gdf, values, spatial_weights, unique_id, rng=None):
         try:
-            from pysal.explore.inequality.theil import Theil
+            from inequality.theil import Theil
         except ImportError:
-            raise ImportError("The 'inequality' or 'pysal' package is required.")
-
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-
-    print("Calculating Theil index...")
-
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
-
-    data = data.set_index(unique_id)
+            try:
+                from pysal.explore.inequality.theil import Theil
+            except ImportError:
+                raise ImportError("The 'inequality' or 'pysal' package is required.")
+
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.rng = rng
+
+        data = gdf.copy()
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
+
+        data = data.set_index(unique_id)
+
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
+
+            values_list = data.loc[neighbours][values]
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = [index]
-        values_list = data.loc[neighbours][values]
-
-        if rng:
-            from momepy import limit_range
+            if rng:
+                from momepy import limit_range
 
-            values_list = limit_range(values_list.tolist(), rng=rng)
-        if set(values_list) == {0}:
-            results_list.append(0)
-        else:
+                values_list = limit_range(values_list, rng=rng)
             results_list.append(Theil(values_list).T)
 
-    series = pd.Series(results_list, index=gdf.index)
+        self.t = pd.Series(results_list, index=gdf.index)
 
-    print("Theil index calculated.")
-    return series
 
-
-def simpson(
-    gdf,
-    values,
-    spatial_weights,
-    unique_id,
-    binning="HeadTailBreaks",
-    **classification_kwds
-):
+class Simpson:
     """
     Calculates the Simpson\'s diversity index of values within neighbours defined in `spatial_weights`.
 
@@ -203,10 +216,22 @@ def simpson(
         For details see mapclassify documentation:
         https://pysal.org/mapclassify
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    s : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    binning : str
+        binning
+    classification_kwds : dict
+        classification_kwds
 
     References
     ----------
@@ -216,13 +241,67 @@ def simpson(
     Examples
     --------
     >>> sw = momepy.sw_high(k=3, gdf=tessellation_df, ids='uID')
-    >>> tessellation_df['area_Simpson'] = mm.simpson(tessellation_df, 'area', sw, 'uID')
-    Calculating Simpson's diversity index...
+    >>> tessellation_df['area_Simpson'] = mm.Simpson(tessellation_df, 'area', sw, 'uID').s
     100%|██████████| 144/144 [00:00<00:00, 455.83it/s]
-    Simpson's diversity index calculated.
     """
 
-    def _simpson_di(data):
+    def __init__(
+        self,
+        gdf,
+        values,
+        spatial_weights,
+        unique_id,
+        binning="HeadTailBreaks",
+        **classification_kwds
+    ):
+        try:
+            import mapclassify.classifiers as classifiers
+        except ImportError:
+            try:
+                import pysal.viz.mapclassify.classifiers as classifiers
+            except ImportError:
+                raise ImportError("The 'mapclassify' or 'pysal' package is required")
+
+        schemes = {}
+        for classifier in classifiers.CLASSIFIERS:
+            schemes[classifier.lower()] = getattr(classifiers, classifier)
+        binning = binning.lower()
+        if binning not in schemes:
+            raise ValueError(
+                "Invalid binning. Binning must be in the" " set: %r" % schemes.keys()
+            )
+
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.binning = binning
+        self.classification_kwds = classification_kwds
+
+        data = gdf.copy()
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
+
+        bins = schemes[binning](data[values], **classification_kwds).bins
+        data = data.set_index(unique_id)
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = [index]
+            values_list = data.loc[neighbours][values]
+
+            sample_bins = classifiers.UserDefined(values_list, bins)
+            counts = dict(zip(bins, sample_bins.counts))
+            results_list.append(self._simpson_di(counts))
+
+        self.s = pd.Series(results_list, index=gdf.index)
+
+    def _simpson_di(self, data):
 
         """ Given a hash { 'species': count } , returns the Simpson Diversity Index
 
@@ -242,55 +321,8 @@ def simpson(
 
         return sum(p(n, N) ** 2 for n in data.values() if n != 0)
 
-    try:
-        import mapclassify.classifiers as classifiers
-    except ImportError:
-        try:
-            import pysal.viz.mapclassify.classifiers as classifiers
-        except ImportError:
-            raise ImportError("The 'mapclassify' or 'pysal' package is required")
-
-    schemes = {}
-    for classifier in classifiers.CLASSIFIERS:
-        schemes[classifier.lower()] = getattr(classifiers, classifier)
-    binning = binning.lower()
-    if binning not in schemes:
-        raise ValueError(
-            "Invalid binning. Binning must be in the" " set: %r" % schemes.keys()
-        )
-
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-    print("Calculating Simpson's diversity index...")
-
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
-
-    bins = schemes[binning](data[values], **classification_kwds).bins
-    data = data.set_index(unique_id)
-
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = [index]
-        values_list = data.loc[neighbours][values]
-
-        sample_bins = classifiers.UserDefined(values_list, bins)
-        counts = dict(zip(bins, sample_bins.counts))
-        results_list.append(_simpson_di(counts))
-
-    series = pd.Series(results_list, index=gdf.index)
-
-    print("Simpson's diversity index calculated.")
-    return series
 
-
-def gini(gdf, values, spatial_weights, unique_id, rng=None):
+class Gini:
     """
     Calculates the Gini index of values within neighbours defined in `spatial_weights`.
 
@@ -310,58 +342,71 @@ def gini(gdf, values, spatial_weights, unique_id, rng=None):
         Percentiles over which to compute the range. Each must be
         between 0 and 100, inclusive. The order of the elements is not important.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    g : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    rng : tuple
+        range
 
     Examples
     --------
     >>> sw = momepy.sw_high(k=3, gdf=tessellation_df, ids='uID')
-    >>> tessellation_df['area_Gini'] = mm.gini(tessellation_df, 'area', sw, 'uID')
-    Calculating Gini index...
+    >>> tessellation_df['area_Gini'] = mm.Gini(tessellation_df, 'area', sw, 'uID').g
     100%|██████████| 144/144 [00:00<00:00, 597.37it/s]
-    Gini index calculated.
     """
-    try:
-        from inequality.gini import Gini
-    except ImportError:
+
+    def __init__(self, gdf, values, spatial_weights, unique_id, rng=None):
         try:
-            from pysal.explore.inequality.gini import Gini
+            from inequality.gini import Gini
         except ImportError:
-            raise ImportError("The 'inequality' or 'pysal' package is required.")
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
-
-    if data[values].min() < 0:
-        raise ValueError(
-            "Values contain negative numbers. Normalise data before"
-            "using momepy.gini."
-        )
-    print("Calculating Gini index...")
-
-    data = data.set_index(unique_id)
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-
-            values_list = data.loc[neighbours][values].values
-
-            if rng:
-                from momepy import limit_range
-
-                values_list = np.array(limit_range(values_list, rng=rng))
-            results_list.append(Gini(values_list).g)
-        else:
-            results_list.append(0)
-    series = pd.Series(results_list, index=gdf.index)
-
-    print("Gini index calculated.")
-    return series
+            try:
+                from pysal.explore.inequality.gini import Gini
+            except ImportError:
+                raise ImportError("The 'inequality' or 'pysal' package is required.")
+
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.rng = rng
+
+        data = gdf.copy()
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
+
+        if self.values.min() < 0:
+            raise ValueError(
+                "Values contain negative numbers. Normalise data before"
+                "using momepy.Gini."
+            )
+
+        data = data.set_index(unique_id)
+
+        results_list = []
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+
+                values_list = data.loc[neighbours][values].values
+
+                if rng:
+                    from momepy import limit_range
+
+                    values_list = np.array(limit_range(values_list, rng=rng))
+                results_list.append(Gini(values_list).g)
+            else:
+                results_list.append(0)
+
+        self.g = pd.Series(results_list, index=gdf.index)
diff --git a/momepy/intensity.py b/momepy/intensity.py
index e44699d..ccab0db 100644
--- a/momepy/intensity.py
+++ b/momepy/intensity.py
@@ -10,15 +10,7 @@ import numpy as np
 import collections
 
 
-def object_area_ratio(
-    left,
-    right,
-    left_areas,
-    right_areas,
-    unique_id=None,
-    left_unique_id=None,
-    right_unique_id=None,
-):
+class AreaRatio:
     """
     Calculate covered area ratio or floor area ratio of objects.
 
@@ -46,10 +38,22 @@ def object_area_ratio(
     right_unique_id : str, list, np.array, pd.Series (default None)
         the name of the left dataframe column, np.array, or pd.Series where is stored shared unique ID
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ar : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    left_areas :  Series
+        Series containing used left areas
+    right_areas :  Series
+        Series containing used right areas
+    left_unique_id : Series
+        Series containing used left ID
+    right_unique_id : Series
+        Series containing used right ID
 
     References
     ---------
@@ -58,51 +62,60 @@ def object_area_ratio(
 
     Examples
     --------
-    >>> tessellation_df['CAR'] = mm.object_area_ratio(tessellation_df, buildings_df, 'area', 'area', 'uID')
-    Calculating object area ratio...
-    Merging DataFrames...
-    Calculating OAR...
-    Object area ratio calculated.
+    >>> tessellation_df['CAR'] = mm.AreaRatio(tessellation_df, buildings_df, 'area', 'area', 'uID').ar
     """
-    left = left.copy()
-    right = right.copy()
-
-    if unique_id:
-        left_unique_id = unique_id
-        right_unique_id = unique_id
-    else:
-        if left_unique_id is None or right_unique_id is None:
-            raise ValueError(
-                "Unique ID not correctly set. Use either network_id or both"
-                "left_unique_id and right_unique_id."
-            )
-
-    print("Calculating object area ratio...")
-    if not isinstance(left_areas, str):
-        left["mm_a"] = left_areas
-        left_areas = "mm_a"
-    if not isinstance(right_areas, str):
-        right["mm_a"] = right_areas
-        right_areas = "mm_a"
-
-    print("Merging DataFrames...")
-    look_for = right[[right_unique_id, right_areas]]  # keeping only necessary columns
-    look_for.rename(index=str, columns={right_areas: "lf_area"}, inplace=True)
-    objects_merged = left[[left_unique_id, left_areas]].merge(
-        look_for, left_on=left_unique_id, right_on=right_unique_id
-    )
-
-    print("Calculating OAR...")
-
-    series = objects_merged["lf_area"] / objects_merged[left_areas]
-
-    print("Object area ratio calculated.")
-    return series
-
-
-def elements_count(left, right, left_id, right_id, weighted=False):
+
+    def __init__(
+        self,
+        left,
+        right,
+        left_areas,
+        right_areas,
+        unique_id=None,
+        left_unique_id=None,
+        right_unique_id=None,
+    ):
+        self.left = left
+        self.right = right
+
+        left = left.copy()
+        right = right.copy()
+
+        if unique_id:
+            left_unique_id = unique_id
+            right_unique_id = unique_id
+        else:
+            if left_unique_id is None or right_unique_id is None:
+                raise ValueError(
+                    "Unique ID not correctly set. Use either network_id or both"
+                    "left_unique_id and right_unique_id."
+                )
+        self.left_unique_id = left_unique_id
+        self.right_unique_id = right_unique_id
+
+        if not isinstance(left_areas, str):
+            left["mm_a"] = left_areas
+            left_areas = "mm_a"
+        self.left_areas = left[left_areas]
+        if not isinstance(right_areas, str):
+            right["mm_a"] = right_areas
+            right_areas = "mm_a"
+        self.right_areas = right[right_areas]
+
+        look_for = right[
+            [right_unique_id, right_areas]
+        ]  # keeping only necessary columns
+        look_for.rename(index=str, columns={right_areas: "lf_area"}, inplace=True)
+        objects_merged = left[[left_unique_id, left_areas]].merge(
+            look_for, left_on=left_unique_id, right_on=right_unique_id
+        )
+
+        self.ar = objects_merged["lf_area"] / objects_merged[left_areas]
+
+
+class Count:
     """
-    Calculate the number of elements within aggregated structure.
+    Calculate the number of elements within an aggregated structure.
 
     Aggregated structure can be typically block, street segment or street node (their snapepd objects). Right gdf has to have
     unique id of aggregated structure assigned before hand (e.g. using :py:func:`momepy.get_network_id`).
@@ -125,10 +138,20 @@ def elements_count(left, right, left_id, right_id, weighted=False):
     weighted : bool (default False)
         if weighted=True, count will be divided by the area or length
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    c : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    left_id : Series
+        Series containing used left ID
+    right_id : Series
+        Series containing used right ID
+    weighted : bool
+        used weighted value
 
     References
     ----------
@@ -140,25 +163,33 @@ def elements_count(left, right, left_id, right_id, weighted=False):
 
     Examples
     --------
-    >>> blocks_df['buildings_count'] = mm.elements_count(blocks_df, buildings_df, 'bID', 'bID', weighted=True)
+    >>> blocks_df['buildings_count'] = mm.Count(blocks_df, buildings_df, 'bID', 'bID', weighted=True).c
     """
-    count = collections.Counter(right[right_id])
-    df = pd.DataFrame.from_dict(count, orient="index", columns=["mm_count"])
-    joined = left[[left_id, "geometry"]].join(df["mm_count"], on=left_id)
-    joined.loc[joined["mm_count"].isna(), "mm_count"] = 0
-
-    if weighted:
-        if left.geometry[0].type in ["Polygon", "MultiPolygon"]:
-            joined["mm_count"] = joined["mm_count"] / left.geometry.area
-        elif left.geometry[0].type in ["LineString", "MultiLineString"]:
-            joined["mm_count"] = joined["mm_count"] / left.geometry.length
-        else:
-            raise TypeError("Geometry type does not support weighting.")
 
-    return joined["mm_count"]
+    def __init__(self, left, right, left_id, right_id, weighted=False):
+        self.left = left
+        self.right = right
+        self.left_id = left[left_id]
+        self.right_id = right[right_id]
+        self.weighted = weighted
+
+        count = collections.Counter(right[right_id])
+        df = pd.DataFrame.from_dict(count, orient="index", columns=["mm_count"])
+        joined = left[[left_id, "geometry"]].join(df["mm_count"], on=left_id)
+        joined.loc[joined["mm_count"].isna(), "mm_count"] = 0
+
+        if weighted:
+            if left.geometry[0].type in ["Polygon", "MultiPolygon"]:
+                joined["mm_count"] = joined["mm_count"] / left.geometry.area
+            elif left.geometry[0].type in ["LineString", "MultiLineString"]:
+                joined["mm_count"] = joined["mm_count"] / left.geometry.length
+            else:
+                raise TypeError("Geometry type does not support weighting.")
+
+        self.c = joined["mm_count"]
 
 
-def courtyards(gdf, block_id, spatial_weights=None):
+class Courtyards:
     """
     Calculate the number of courtyards within the joined structure.
 
@@ -172,10 +203,16 @@ def courtyards(gdf, block_id, spatial_weights=None):
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on objects. It is to denote adjacent buildings (note: based on index).
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    c : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    block_id : Series
+        Series containing used block ID
+    sw : libpysal.weights
+        spatial weights matrix
 
     References
     ---------
@@ -184,62 +221,57 @@ def courtyards(gdf, block_id, spatial_weights=None):
 
     Examples
     --------
-    >>> buildings_df['courtyards'] = mm.courtyards(buildings_df, 'bID')
-    Calculating courtyards...
+    >>> buildings_df['courtyards'] = mm.Courtyards(buildings_df, 'bID').c
     Calculating spatial weights...
-    Spatial weights ready...
-    Courtyards calculated.
-
-    Notes
-    -----
-    Script is not well optimised, so it is currently slow.
     """
-    # define empty list for results
-    results_list = []
-
-    print("Calculating courtyards...")
-    if not isinstance(block_id, str):
-        gdf["mm_bid"] = block_id
-        block_id = "mm_bid"
-
-    # if weights matrix is not passed, generate it from objects
-    if spatial_weights is None:
-        print("Calculating spatial weights...")
-        from libpysal.weights import Queen
-
-        spatial_weights = Queen.from_dataframe(gdf, silence_warnings=True)
-        print("Spatial weights ready...")
-
-    # dict to store nr of courtyards for each uID
-    courtyards = {}
-    components = pd.Series(spatial_weights.component_labels, index=gdf.index)
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        # if the id is already present in courtyards, continue (avoid repetition)
-        if index in courtyards:
-            continue
-        else:
-            comp = spatial_weights.component_labels[index]
-            to_join = components[components == comp].index
-            joined = gdf.loc[to_join]
-            dissolved = joined.geometry.buffer(
-                0.01
-            ).unary_union  # buffer to avoid multipolygons where buildings touch by corners only
-            try:
-                interiors = len(list(dissolved.interiors))
-            except (ValueError):
-                print("Something unexpected happened.")
-            for b in to_join:
-                courtyards[b] = interiors  # fill dict with values
-    # copy values from dict to gdf
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        results_list.append(courtyards[index])
-
-    series = pd.Series(results_list, index=gdf.index)
-    print("Courtyards calculated.")
-    return series
-
-
-def blocks_count(gdf, block_id, spatial_weights, unique_id, weighted=True):
+
+    def __init__(self, gdf, block_id, spatial_weights=None):
+        self.gdf = gdf
+
+        results_list = []
+        gdf = gdf.copy()
+
+        if not isinstance(block_id, str):
+            gdf["mm_bid"] = block_id
+            block_id = "mm_bid"
+
+        self.block_id = gdf[block_id]
+        # if weights matrix is not passed, generate it from objects
+        if spatial_weights is None:
+            print("Calculating spatial weights...")
+            from libpysal.weights import Queen
+
+            spatial_weights = Queen.from_dataframe(gdf, silence_warnings=True)
+
+        self.sw = spatial_weights
+        # dict to store nr of courtyards for each uID
+        courtyards = {}
+        components = pd.Series(spatial_weights.component_labels, index=gdf.index)
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            # if the id is already present in courtyards, continue (avoid repetition)
+            if index in courtyards:
+                continue
+            else:
+                comp = spatial_weights.component_labels[index]
+                to_join = components[components == comp].index
+                joined = gdf.loc[to_join]
+                dissolved = joined.geometry.buffer(
+                    0.01
+                ).unary_union  # buffer to avoid multipolygons where buildings touch by corners only
+                try:
+                    interiors = len(list(dissolved.interiors))
+                except (ValueError):
+                    print("Something unexpected happened.")
+                for b in to_join:
+                    courtyards[b] = interiors  # fill dict with values
+        # copy values from dict to gdf
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            results_list.append(courtyards[index])
+
+        self.c = pd.Series(results_list, index=gdf.index)
+
+
+class BlocksCount:
     """
     Calculates the weighted number of blocks
 
@@ -261,10 +293,20 @@ def blocks_count(gdf, block_id, spatial_weights, unique_id, weighted=True):
     weigted : bool, default True
         return value weighted by the analysed area (True) or pure count (False)
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    bc : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    block_id : Series
+        Series containing used block ID
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    weighted : bool
+        used weighted value
 
     References
     ----------
@@ -274,47 +316,46 @@ def blocks_count(gdf, block_id, spatial_weights, unique_id, weighted=True):
     Examples
     --------
     >>> sw4 = mm.sw_high(k=4, gdf='tessellation_df', ids='uID')
-    >>> tessellation_df['blocks_within_4'] = mm.blocks(tessellation_df, 'bID', sw4, 'uID')
-    Calculating blocks...
-    Blocks calculated.
+    >>> tessellation_df['blocks_within_4'] = mm.BlocksCount(tessellation_df, 'bID', sw4, 'uID').bc
     """
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-    if not isinstance(block_id, str):
-        data["mm_bid"] = block_id
-        block_id = "mm_bid"
-
-    data = data.set_index(unique_id)
-
-    print("Calculating blocks...")
 
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = row[unique_id]
-        vicinity = data.loc[neighbours]
-
-        if weighted is True:
-            results_list.append(
-                len(set(list(vicinity[block_id]))) / sum(vicinity.geometry.area)
-            )
-        elif weighted is False:
-            results_list.append(len(set(list(vicinity[block_id]))))
-        else:
-            raise ValueError("Attribute 'weighted' needs to be True or False.")
+    def __init__(self, gdf, block_id, spatial_weights, unique_id, weighted=True):
+
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+        self.weighted = weighted
+
+        # define empty list for results
+        results_list = []
+        data = gdf.copy()
+        if not isinstance(block_id, str):
+            data["mm_bid"] = block_id
+            block_id = "mm_bid"
+        self.block_id = data[block_id]
+        data = data.set_index(unique_id)
+
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = row[unique_id]
+            vicinity = data.loc[neighbours]
 
-    series = pd.Series(results_list, index=gdf.index)
+            if weighted is True:
+                results_list.append(
+                    len(set(list(vicinity[block_id]))) / sum(vicinity.geometry.area)
+                )
+            elif weighted is False:
+                results_list.append(len(set(list(vicinity[block_id]))))
+            else:
+                raise ValueError("Attribute 'weighted' needs to be True or False.")
 
-    print("Blocks calculated.")
-    return series
+        self.bc = pd.Series(results_list, index=gdf.index)
 
 
-def reached(
-    left, right, left_id, right_id, spatial_weights=None, mode="count", values=None
-):
+class Reached:
     """
     Calculates the number of objects reached within topological steps on street network
 
@@ -349,90 +390,105 @@ def reached(
     values : str (default None)
         the name of the objects dataframe column with values used for calculations
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    r : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    left_id : Series
+        Series containing used left ID
+    right_id : Series
+        Series containing used right ID
+    mode : str
+        mode of calculation
+    sw : libpysal.weights
+        spatial weights matrix (if set)
 
     Examples
     --------
-    >>> streets_df['reached_buildings'] = mm.reached(streets_df, buildings_df, 'uID')
-    Calculating reached count...
-    Reached count calculated.
+    >>> streets_df['reached_buildings'] = mm.Reached(streets_df, buildings_df, 'uID').r
 
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating reached {}...".format(mode))
-
-    if not isinstance(right_id, str):
-        right = right.copy()
-        right["mm_id"] = right_id
-        right_id = "mm_id"
-
-    if not isinstance(left_id, str):
-        left = left.copy()
-        left["mm_lid"] = left_id
-        left_id = "mm_lid"
-
-    if mode == "count":
-        count = collections.Counter(right[right_id])
-
-    # iterating over rows one by one
-    for index, row in tqdm(left.iterrows(), total=left.shape[0]):
-        if spatial_weights is None:
-            ids = [row[left_id]]
-        else:
-            neighbours = list(spatial_weights.neighbors[index])
-            neighbours.append(index)
-            ids = left.iloc[neighbours][left_id]
+    def __init__(
+        self,
+        left,
+        right,
+        left_id,
+        right_id,
+        spatial_weights=None,
+        mode="count",
+        values=None,
+    ):
+        self.left = left
+        self.right = right
+        self.sw = spatial_weights
+        self.mode = mode
+
+        # define empty list for results
+        results_list = []
+
+        if not isinstance(right_id, str):
+            right = right.copy()
+            right["mm_id"] = right_id
+            right_id = "mm_id"
+        self.right_id = right[right_id]
+        if not isinstance(left_id, str):
+            left = left.copy()
+            left["mm_lid"] = left_id
+            left_id = "mm_lid"
+        self.left_id = left[left_id]
         if mode == "count":
-            counts = []
-            for nid in ids:
-                counts.append(count[nid])
-            results_list.append(sum(counts))
-        elif mode == "sum":
-            if values:
-                results_list.append(sum(right.loc[right[right_id].isin(ids)][values]))
-            else:
-                results_list.append(
-                    sum(right.loc[right[right_id].isin(ids)].geometry.area)
-                )
-        elif mode == "mean":
-            if values:
-                results_list.append(
-                    np.nanmean(right.loc[right[right_id].isin(ids)][values])
-                )
-            else:
-                results_list.append(
-                    np.nanmean(right.loc[right[right_id].isin(ids)].geometry.area)
-                )
-        elif mode == "std":
-            if values:
-                results_list.append(
-                    np.nanstd(right.loc[right[right_id].isin(ids)][values])
-                )
-            else:
-                results_list.append(
-                    np.nanstd(right.loc[right[right_id].isin(ids)].geometry.area)
-                )
-
-    series = pd.Series(results_list, index=left.index)
-
-    print("Reached {} calculated.".format(mode))
-    return series
-
+            count = collections.Counter(right[right_id])
 
-def node_density(
-    left,
-    right,
-    spatial_weights,
-    weighted=False,
-    node_degree=None,
-    node_start="node_start",
-    node_end="node_end",
-):
+        # iterating over rows one by one
+        for index, row in tqdm(left.iterrows(), total=left.shape[0]):
+            if spatial_weights is None:
+                ids = [row[left_id]]
+            else:
+                neighbours = list(spatial_weights.neighbors[index])
+                neighbours.append(index)
+                ids = left.iloc[neighbours][left_id]
+            if mode == "count":
+                counts = []
+                for nid in ids:
+                    counts.append(count[nid])
+                results_list.append(sum(counts))
+            elif mode == "sum":
+                if values:
+                    results_list.append(
+                        sum(right.loc[right[right_id].isin(ids)][values])
+                    )
+                else:
+                    results_list.append(
+                        sum(right.loc[right[right_id].isin(ids)].geometry.area)
+                    )
+            elif mode == "mean":
+                if values:
+                    results_list.append(
+                        np.nanmean(right.loc[right[right_id].isin(ids)][values])
+                    )
+                else:
+                    results_list.append(
+                        np.nanmean(right.loc[right[right_id].isin(ids)].geometry.area)
+                    )
+            elif mode == "std":
+                if values:
+                    results_list.append(
+                        np.nanstd(right.loc[right[right_id].isin(ids)][values])
+                    )
+                else:
+                    results_list.append(
+                        np.nanstd(right.loc[right[right_id].isin(ids)].geometry.area)
+                    )
+
+        self.r = pd.Series(results_list, index=left.index)
+
+
+class NodeDensity:
     """
     Calculate the density of nodes within topological steps on street network defined in spatial_weights.
 
@@ -448,21 +504,36 @@ def node_density(
         GeoDataFrame containing nodes of street network
     right : GeoDataFrame
         GeoDataFrame containing edges of street network
-    spatial_weights : libpysal.weights, optional
+    spatial_weights : libpysal.weights
         spatial weights matrix capturing relationship between nodes within set topological distance
-    weighted : bool
+    weighted : bool (default False)
         if True density will take into account node degree as k-1
-    node_degree : str
-        name of the column of left gdf containing node degree
+    node_degree : str (optional)
+        name of the column of left gdf containing node degree. Used if weighted=True
     node_start : str (default 'node_start')
         name of the column of right gdf containing id of starting node
     node_end : str (default 'node_end')
         name of the column of right gdf containing id of ending node
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    nd : Series
+        Series containing resulting values
+    left : GeoDataFrame
+        original left GeoDataFrame
+    right : GeoDataFrame
+        original right GeoDataFrame
+    node_start : Series
+        Series containing used ids of starting node
+    node_end : Series
+        Series containing used ids of ending node
+    sw : libpysal.weights
+        spatial weights matrix
+    weighted : bool
+        used weighted value
+    node_degree : Series
+        Series containing used node degree values
+
 
     References
     ---------
@@ -471,43 +542,56 @@ def node_density(
 
     Examples
     --------
-    >>> nodes['density'] = mm.node_density(nodes, edges, sw)
-    Calculating node density...
-    Node density calculated.
+    >>> nodes['density'] = mm.NodeDensity(nodes, edges, sw).nd
 
     """
-    # define empty list for results
-    results_list = []
 
-    print("Calculating node density...")
+    def __init__(
+        self,
+        left,
+        right,
+        spatial_weights,
+        weighted=False,
+        node_degree=None,
+        node_start="node_start",
+        node_end="node_end",
+    ):
+        self.left = left
+        self.right = right
+        self.sw = spatial_weights
+        self.weighted = weighted
+        if weighted:
+            self.node_degree = left[node_degree]
+        self.node_start = right[node_start]
+        self.node_end = right[node_end]
+        # define empty list for results
+        results_list = []
 
-    # iterating over rows one by one
-    for index, row in tqdm(left.iterrows(), total=left.shape[0]):
+        # iterating over rows one by one
+        for index, row in tqdm(left.iterrows(), total=left.shape[0]):
 
-        neighbours = list(spatial_weights.neighbors[index])
-        neighbours.append(index)
-        if weighted:
-            neighbour_nodes = left.iloc[neighbours]
-            number_nodes = sum(neighbour_nodes[node_degree] - 1)
-        else:
-            number_nodes = len(neighbours)
+            neighbours = list(spatial_weights.neighbors[index])
+            neighbours.append(index)
+            if weighted:
+                neighbour_nodes = left.iloc[neighbours]
+                number_nodes = sum(neighbour_nodes[node_degree] - 1)
+            else:
+                number_nodes = len(neighbours)
 
-        edg = right.loc[right["node_start"].isin(neighbours)].loc[
-            right["node_end"].isin(neighbours)
-        ]
-        length = sum(edg.geometry.length)
+            edg = right.loc[right["node_start"].isin(neighbours)].loc[
+                right["node_end"].isin(neighbours)
+            ]
+            length = sum(edg.geometry.length)
 
-        if length > 0:
-            results_list.append(number_nodes / length)
-        else:
-            results_list.append(0)
+            if length > 0:
+                results_list.append(number_nodes / length)
+            else:
+                results_list.append(0)
 
-    series = pd.Series(results_list, index=left.index)
-    print("Node density calculated.")
-    return series
+        self.nd = pd.Series(results_list, index=left.index)
 
 
-def density(gdf, values, spatial_weights, unique_id, areas=None):
+class Density:
     """
     Calculate the gross density
 
@@ -524,14 +608,24 @@ def density(gdf, values, spatial_weights, unique_id, areas=None):
         spatial weights matrix
     unique_id : str
         name of the column with unique id used as spatial_weights index
-    areas :  str, list, np.array, pd.Series
+    areas :  str, list, np.array, pd.Series (optional)
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If None,
         gdf.geometry.area will be used.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    d : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    values : Series
+        Series containing used values
+    sw : libpysal.weights
+        spatial weights matrix
+    id : Series
+        Series containing used unique ID
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -540,40 +634,44 @@ def density(gdf, values, spatial_weights, unique_id, areas=None):
 
     Examples
     --------
-    >>> tessellation_df['floor_area_dens'] = mm.density(tessellation_df, 'floor_area', sw, 'uID')
-    Calculating gross density...
-    Gross density calculated.
+    >>> tessellation_df['floor_area_dens'] = mm.Density(tessellation_df, 'floor_area', sw, 'uID').d
     """
-    # define empty list for results
-    results_list = []
-    data = gdf.copy()
-
-    print("Calculating gross density...")
-    if values is not None:
-        if not isinstance(values, str):
-            data["mm_v"] = values
-            values = "mm_v"
-    if areas is not None:
-        if not isinstance(areas, str):
-            data["mm_a"] = areas
-            areas = "mm_a"
-    data = data.set_index(unique_id)
-    # iterating over rows one by one
-    for index, row in tqdm(data.iterrows(), total=data.shape[0]):
-        neighbours = spatial_weights.neighbors[index].copy()
-        if neighbours:
-            neighbours.append(index)
-        else:
-            neighbours = index
-        subset = data.loc[neighbours]
-        values_list = subset[values]
+
+    def __init__(self, gdf, values, spatial_weights, unique_id, areas=None):
+        self.gdf = gdf
+        self.sw = spatial_weights
+        self.id = gdf[unique_id]
+
+        # define empty list for results
+        results_list = []
+        data = gdf.copy()
+
+        if values is not None:
+            if not isinstance(values, str):
+                data["mm_v"] = values
+                values = "mm_v"
+        self.values = data[values]
         if areas is not None:
-            areas_list = subset[areas]
+            if not isinstance(areas, str):
+                data["mm_a"] = areas
+                areas = "mm_a"
         else:
-            areas_list = subset.geometry.area
+            data["mm_a"] = data.geometry.area
+            areas = "mm_a"
+        self.areas = data[areas]
+
+        data = data.set_index(unique_id)
+        # iterating over rows one by one
+        for index, row in tqdm(data.iterrows(), total=data.shape[0]):
+            neighbours = spatial_weights.neighbors[index].copy()
+            if neighbours:
+                neighbours.append(index)
+            else:
+                neighbours = index
+            subset = data.loc[neighbours]
+            values_list = subset[values]
+            areas_list = subset[areas]
 
-        results_list.append(sum(values_list) / sum(areas_list))
+            results_list.append(sum(values_list) / sum(areas_list))
 
-    series = pd.Series(results_list, index=gdf.index)
-    print("Gross density calculated.")
-    return series
+        self.d = pd.Series(results_list, index=gdf.index)
diff --git a/momepy/shape.py b/momepy/shape.py
index ba37926..8afe4b3 100644
--- a/momepy/shape.py
+++ b/momepy/shape.py
@@ -12,7 +12,7 @@ import numpy as np
 from shapely.geometry import Point
 
 
-def form_factor(gdf, volumes, areas=None):
+class FormFactor:
     """
     Calculates form factor of each object in given geoDataFrame.
 
@@ -30,10 +30,16 @@ def form_factor(gdf, volumes, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ff : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    volumes : Series
+        Series containing used volume values
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -43,50 +49,39 @@ def form_factor(gdf, volumes, areas=None):
 
     Examples
     --------
-    >>> buildings_df['formfactor'] = momepy.form_factor(buildings_df, 'volume')
-    Calculating form factor...
-    Form factor calculated.
+    >>> buildings_df['formfactor'] = momepy.FormFactor(buildings_df, 'volume').ff
     >>> buildings_df.formfactor[0]
     1.9385988170288635
 
-    >>> buildings_df['formfactor'] = momepy.form_factor(buildings_df, momepy.volume(buildings_df, 'height'))
-    Calculating volumes...
-    Volumes calculated.
-    Calculating form factor...
-    Form factor calculated.
+    >>> buildings_df['formfactor'] = momepy.FormFactor(buildings_df, momepy.volume(buildings_df, 'height').volume).ff
     >>> buildings_df.formfactor[0]
     1.9385988170288635
 
     """
 
-    print("Calculating form factor...")
-    gdf = gdf.copy()
-    if not isinstance(volumes, str):
-        gdf["mm_v"] = volumes
-        volumes = "mm_v"
-    if areas is None:
-        series = gdf.apply(
-            lambda row: row.geometry.area / (row[volumes] ** (2 / 3))
-            if row[volumes] != 0
-            else 0,
-            axis=1,
-        )
-    else:
+    def __init__(self, gdf, volumes, areas=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        if not isinstance(volumes, str):
+            gdf["mm_v"] = volumes
+            volumes = "mm_v"
+        self.volumes = gdf[volumes]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.ff = gdf.apply(
             lambda row: row[areas] / (row[volumes] ** (2 / 3))
             if row[volumes] != 0
             else 0,
             axis=1,
         )
 
-    print("Form factor calculated.")
-    return series
-
 
-def fractal_dimension(gdf, areas=None, perimeters=None):
+class FractalDimension:
     """
     Calculates fractal dimension of each object in given geoDataFrame.
 
@@ -104,10 +99,16 @@ def fractal_dimension(gdf, areas=None, perimeters=None):
         the name of the dataframe column, np.array, or pd.Series where is stored perimeter value. If set to None, function will calculate perimeters
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    fd : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    perimeters : Series
+        Series containing used perimeter values
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -117,42 +118,37 @@ def fractal_dimension(gdf, areas=None, perimeters=None):
 
     Examples
     --------
-    >>> buildings_df['fractal'] = momepy.fractal_dimension(buildings_df, 'area', 'peri')
-    Calculating fractal dimension...
+    >>> buildings_df['fractal'] = momepy.FractalDimension(buildings_df, 'area', 'peri').fd
     100%|██████████| 144/144 [00:00<00:00, 3928.09it/s]
-    Fractal dimension calculated.
     >>> buildings_df.fractal[0]
     0.5363389283519454
     """
 
-    print("Calculating fractal dimension...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None, perimeters=None):
+        self.gdf = gdf
 
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+        gdf = gdf.copy()
+
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
-    if areas is None:
-        series = gdf.apply(
-            lambda row: math.log(row[perimeters] / 4) / math.log(row.geometry.area),
-            axis=1,
-        )
-    else:
+        else:
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+        self.perimeters = gdf[perimeters]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+
+        self.fd = gdf.apply(
             lambda row: math.log(row[perimeters] / 4) / math.log(row[areas]), axis=1
         )
 
-    print("Fractal dimension calculated.")
-    return series
-
 
-def volume_facade_ratio(gdf, heights, volumes=None, perimeters=None):
+class VolumeFacadeRatio:
     """
     Calculates volume/facade ratio of each object in given geoDataFrame.
 
@@ -171,10 +167,16 @@ def volume_facade_ratio(gdf, heights, volumes=None, perimeters=None):
         the name of the dataframe column, np.array, or pd.Series where is stored perimeter value
 
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    vfr : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    perimeters : Series
+        Series containing used perimeter values
+    volumes : Series
+        Series containing used volume values
 
     References
     ----------
@@ -184,35 +186,34 @@ def volume_facade_ratio(gdf, heights, volumes=None, perimeters=None):
 
     Examples
     -----
-    >>> buildings_df['vfr'] = momepy.volume_facade_ratio(buildings_df, 'height')
-    Calculating volume/facade ratio...
-    Volume/facade ratio calculated.
+    >>> buildings_df['vfr'] = momepy.VolumeFacadeRatio(buildings_df, 'height').vfr
     >>> buildings_df.vfr[0]
     5.310715735236504
     """
 
-    print("Calculating volume/facade ratio...")
-    gdf = gdf.copy()
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+    def __init__(self, gdf, heights, volumes=None, perimeters=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
+        else:
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+        self.perimeters = gdf[perimeters]
 
-    if volumes is None:
-        gdf["mm_v"] = gdf.geometry.area * gdf[heights]
-        volumes = "mm_v"
-    else:
-        if not isinstance(volumes, str):
-            gdf["mm_v"] = volumes
+        if volumes is None:
+            gdf["mm_v"] = gdf.geometry.area * gdf[heights]
             volumes = "mm_v"
+        else:
+            if not isinstance(volumes, str):
+                gdf["mm_v"] = volumes
+                volumes = "mm_v"
+        self.volumes = gdf[volumes]
 
-    series = gdf[volumes] / (gdf[perimeters] * gdf[heights])
-
-    print("Volume/facade ratio calculated.")
-    return series
+        self.vfr = gdf[volumes] / (gdf[perimeters] * gdf[heights])
 
 
 # Smallest enclosing circle - Library (Python)
@@ -384,7 +385,7 @@ def _circle_area(points):
     return math.pi * circ[2] ** 2
 
 
-def circular_compactness(gdf, areas=None):
+class CircularCompactness:
     """
     Calculates compactness index of each object in given geoDataFrame.
 
@@ -399,10 +400,14 @@ def circular_compactness(gdf, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    cc : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
 
     References
     ----------
@@ -411,38 +416,31 @@ def circular_compactness(gdf, areas=None):
 
     Examples
     --------
-    >>> buildings_df['circ_comp'] = momepy.circular_compactness(buildings_df, 'area')
-    Calculating compactness index...
+    >>> buildings_df['circ_comp'] = momepy.CircularCompactness(buildings_df, 'area').cc
     100%|██████████| 144/144 [00:00<00:00, 2498.75it/s]
-    Compactness index calculated.
     >>> buildings_df['circ_comp'][0]
     0.572145421828038
     """
 
-    print("Calculating compactness index...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None):
+        self.gdf = gdf
 
-    if areas is None:
-        series = gdf.apply(
-            lambda row: (row.geometry.area)
-            / (_circle_area(list(row["geometry"].convex_hull.exterior.coords))),
-            axis=1,
-        )
-    else:
+        gdf = gdf.copy()
+
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.cc = gdf.apply(
             lambda row: (row[areas])
             / (_circle_area(list(row["geometry"].convex_hull.exterior.coords))),
             axis=1,
         )
 
-    print("Compactness index calculated.")
-    return series
 
-
-def square_compactness(gdf, areas=None, perimeters=None):
+class SquareCompactness:
     """
     Calculates compactness index of each object in given geoDataFrame.
 
@@ -462,10 +460,16 @@ def square_compactness(gdf, areas=None, perimeters=None):
         the name of the dataframe column, np.array, or pd.Series where is stored perimeter value. If set to None, function will calculate perimeters
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    sc : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
+    perimeters : Series
+        Series containing used perimeter values
 
     References
     ----------
@@ -474,42 +478,37 @@ def square_compactness(gdf, areas=None, perimeters=None):
 
     Examples
     --------
-    >>> buildings_df['squ_comp'] = momepy.square_compactness(buildings_df)
-    Calculating compactness index...
-    Compactness index calculated.
+    >>> buildings_df['squ_comp'] = momepy.SquareCompactness(buildings_df).sc
     >>> buildings_df['squ_comp'][0]
     0.6193872538650996
 
     """
 
-    print("Calculating compactness index...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None, perimeters=None):
+        self.gdf = gdf
 
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+        gdf = gdf.copy()
+
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
-    if areas is None:
-        series = gdf.apply(
-            lambda row: ((4 * math.sqrt(row.geometry.area)) / (row[perimeters])) ** 2,
-            axis=1,
-        )
-    else:
+        else:
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+        self.perimeters = gdf[perimeters]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.sc = gdf.apply(
             lambda row: ((4 * math.sqrt(row[areas])) / (row[perimeters])) ** 2, axis=1
         )
 
-    print("Compactness index calculated.")
-    return series
 
-
-def convexeity(gdf, areas=None):
+class Convexeity:
     """
     Calculates convexeity index of each object in given geoDataFrame.
 
@@ -524,10 +523,14 @@ def convexeity(gdf, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    c : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -536,29 +539,26 @@ def convexeity(gdf, areas=None):
 
     Examples
     --------
-    >>> buildings_df['convexeity'] = momepy.convexeity(buildings_df)
-    Calculating convexeity...
-    Convexeity calculated.
+    >>> buildings_df['convexeity'] = momepy.Convexeity(buildings_df).c
     >>> buildings_df.convexeity[0]
     0.8151964258521672
     """
 
-    print("Calculating convexeity...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None):
+        self.gdf = gdf
+
+        gdf = gdf.copy()
 
-    if areas is None:
-        series = gdf.geometry.area / gdf.geometry.convex_hull.area
-    else:
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf[areas] / gdf.geometry.convex_hull.area
+        self.areas = gdf[areas]
+        self.c = gdf[areas] / gdf.geometry.convex_hull.area
 
-    print("Convexeity calculated.")
-    return series
 
-
-def courtyard_index(gdf, courtyard_areas, areas=None):
+class CourtyardIndex:
     """
     Calculates courtyard index of each object in given geoDataFrame.
 
@@ -576,10 +576,16 @@ def courtyard_index(gdf, courtyard_areas, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    ci : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    courtyard_areas : Series
+        Series containing used courtyard areas values
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -589,40 +595,33 @@ def courtyard_index(gdf, courtyard_areas, areas=None):
 
     Examples
     --------
-    >>> buildings_df['courtyard_index'] = momepy.courtyard_index(buildings, 'courtyard_area', 'area')
-    Calculating courtyard index...
-    Courtyard index calculated.
+    >>> buildings_df['courtyard_index'] = momepy.CourtyardIndex(buildings, 'courtyard_area', 'area').ci
     >>> buildings_df.courtyard_index[80]
     0.16605915738643523
 
-    >>> buildings_df['courtyard_index2'] = momepy.courtyard_index(buildings_df, momepy.courtyard_area(buildings_df))
-    Calculating courtyard areas...
-    Courtyard areas calculated.
-    Calculating courtyard index...
-    Courtyard index calculated.
+    >>> buildings_df['courtyard_index2'] = momepy.CourtyardIndex(buildings_df, momepy.courtyard_area(buildings_df).ca).ci
     >>> buildings_df.courtyard_index2[80]
     0.16605915738643523
     """
 
-    print("Calculating courtyard index...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, courtyard_areas, areas=None):
+        self.gdf = gdf
+        gdf = gdf.copy()
 
-    if not isinstance(courtyard_areas, str):
-        gdf["mm_ca"] = courtyard_areas
-        courtyard_areas = "mm_ca"
-    if areas is None:
-        series = gdf[courtyard_areas] / gdf.geometry.area
-    else:
+        if not isinstance(courtyard_areas, str):
+            gdf["mm_ca"] = courtyard_areas
+            courtyard_areas = "mm_ca"
+        self.courtyard_areas = gdf[courtyard_areas]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf[courtyard_areas] / gdf[areas]
-
-    print("Courtyard index calculated.")
-    return series
+        self.areas = gdf[areas]
+        self.ci = gdf[courtyard_areas] / gdf[areas]
 
 
-def rectangularity(gdf, areas=None):
+class Rectangularity:
     """
     Calculates rectangularity of each object in given geoDataFrame.
 
@@ -637,10 +636,14 @@ def rectangularity(gdf, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    r : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
 
     References
     ----------
@@ -649,37 +652,28 @@ def rectangularity(gdf, areas=None):
 
     Examples
     --------
-    >>> buildings_df['rectangularity'] = momepy.rectangularity(buildings_df, 'area')
-    Calculating rectangularity...
+    >>> buildings_df['rectangularity'] = momepy.Rectangularity(buildings_df, 'area').r
     100%|██████████| 144/144 [00:00<00:00, 866.62it/s]
-    Rectangularity calculated.
     >>> buildings_df.rectangularity[0]
     0.6942676157646379
     """
 
-    print("Calculating rectangularity...")
-    gdf = gdf.copy()
-    if areas is None:
-        series = gdf.apply(
-            lambda row: row.geometry.area
-            / (row.geometry.minimum_rotated_rectangle.area),
-            axis=1,
-        )
-
-    else:
+    def __init__(self, gdf, areas=None):
+        self.gdf = gdf
+        gdf = gdf.copy()
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.r = gdf.apply(
             lambda row: row[areas] / (row.geometry.minimum_rotated_rectangle.area),
             axis=1,
         )
 
-    print("Rectangularity calculated.")
-    return series
-
 
-def shape_index(gdf, longest_axis, areas=None):
+class ShapeIndex:
     """
     Calculates shape index of each object in given geoDataFrame.
 
@@ -696,10 +690,16 @@ def shape_index(gdf, longest_axis, areas=None):
         the name of the dataframe column, np.array, or pd.Series where is stored area value. If set to None, function will calculate areas
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    si : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    longest_axis : Series
+        Series containing used longest axis values
+    areas : Series
+        Series containing used area values
 
     References
     ---------
@@ -707,40 +707,33 @@ def shape_index(gdf, longest_axis, areas=None):
 
     Examples
     --------
-    >>> buildings_df['shape_index'] = momepy.shape_index(buildings_df, longest_axis='long_ax', areas='area')
-    Calculating shape index...
+    >>> buildings_df['shape_index'] = momepy.ShapeIndex(buildings_df, longest_axis='long_ax', areas='area').si
     100%|██████████| 144/144 [00:00<00:00, 5558.33it/s]
-    Shape index calculated.
     >>> buildings_df['shape_index'][0]
     0.7564029493781987
     """
-    print("Calculating shape index...")
-    gdf = gdf.copy()
 
-    if not isinstance(longest_axis, str):
-        gdf["mm_la"] = longest_axis
-        longest_axis = "mm_la"
+    def __init__(self, gdf, longest_axis, areas=None):
+        self.gdf = gdf
+        gdf = gdf.copy()
 
-    if areas is None:
-        series = gdf.apply(
-            lambda row: math.sqrt(row.geometry.area / math.pi)
-            / (0.5 * row[longest_axis]),
-            axis=1,
-        )
-    else:
+        if not isinstance(longest_axis, str):
+            gdf["mm_la"] = longest_axis
+            longest_axis = "mm_la"
+        self.longest_axis = gdf[longest_axis]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.si = gdf.apply(
             lambda row: math.sqrt(row[areas] / math.pi) / (0.5 * row[longest_axis]),
             axis=1,
         )
 
-    print("Shape index calculated.")
-    return series
-
 
-def corners(gdf):
+class Corners:
     """
     Calculates number of corners of each object in given geoDataFrame.
 
@@ -754,83 +747,80 @@ def corners(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    c : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
-    References
-    ---------
-    ?
 
     Examples
     --------
-    >>> buildings_df['corners'] = momepy.corners(buildings_df)
-    Calculating corners...
+    >>> buildings_df['corners'] = momepy.Corners(buildings_df).c
     100%|██████████| 144/144 [00:00<00:00, 1042.15it/s]
-    Corners calculated.
     >>> buildings_df.corners[0]
     24
 
 
     """
-    # define empty list for results
-    results_list = []
-    print("Calculating corners...")
-
-    # calculate angle between points, return true or false if real corner
-    def true_angle(a, b, c):
-        ba = a - b
-        bc = c - b
-
-        cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
-        angle = np.arccos(cosine_angle)
-
-        if np.degrees(angle) <= 170:
-            return True
-        if np.degrees(angle) >= 190:
-            return True
-        return False
-
-    # fill new column with the value of area, iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        corners = 0  # define empty variables
-        points = list(row["geometry"].exterior.coords)  # get points of a shape
-        stop = len(points) - 1  # define where to stop
-        for i in np.arange(
-            len(points)
-        ):  # for every point, calculate angle and add 1 if True angle
-            if i == 0:
-                continue
-            elif i == stop:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[1])
-
-                if true_angle(a, b, c) is True:
-                    corners = corners + 1
-                else:
+
+    def __init__(self, gdf):
+        self.gdf = gdf
+
+        # define empty list for results
+        results_list = []
+
+        # calculate angle between points, return true or false if real corner
+        def _true_angle(a, b, c):
+            ba = a - b
+            bc = c - b
+
+            cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
+            angle = np.arccos(cosine_angle)
+
+            if np.degrees(angle) <= 170:
+                return True
+            if np.degrees(angle) >= 190:
+                return True
+            return False
+
+        # fill new column with the value of area, iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            corners = 0  # define empty variables
+            points = list(row["geometry"].exterior.coords)  # get points of a shape
+            stop = len(points) - 1  # define where to stop
+            for i in np.arange(
+                len(points)
+            ):  # for every point, calculate angle and add 1 if True angle
+                if i == 0:
                     continue
+                elif i == stop:
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[1])
 
-            else:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[i + 1])
+                    if _true_angle(a, b, c) is True:
+                        corners = corners + 1
+                    else:
+                        continue
 
-                if true_angle(a, b, c) is True:
-                    corners = corners + 1
                 else:
-                    continue
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[i + 1])
 
-        results_list.append(corners)
+                    if _true_angle(a, b, c) is True:
+                        corners = corners + 1
+                    else:
+                        continue
 
-    series = pd.Series(results_list, index=gdf.index)
+            results_list.append(corners)
 
-    print("Corners calculated.")
-    return series
+        self.c = pd.Series(results_list, index=gdf.index)
 
 
-def squareness(gdf):
+class Squareness:
     """
     Calculates squareness of each object in given geoDataFrame.
 
@@ -844,10 +834,12 @@ def squareness(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    s : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     References
     ----------
@@ -856,74 +848,71 @@ def squareness(gdf):
 
     Examples
     --------
-    >>> buildings_df['squareness'] = momepy.squareness(buildings_df)
-    Calculating squareness...
+    >>> buildings_df['squareness'] = momepy.Squareness(buildings_df).s
     100%|██████████| 144/144 [00:01<00:00, 129.49it/s]
-    Squareness calculated.
     >>> buildings_df.squareness[0]
     3.7075816043359864
     """
-    # define empty list for results
-    results_list = []
-    print("Calculating squareness...")
-
-    def _angle(a, b, c):
-        ba = a - b
-        bc = c - b
-
-        cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
-        angle = np.degrees(np.arccos(cosine_angle))
-
-        return angle
-
-    # fill new column with the value of area, iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        angles = []
-        points = list(row["geometry"].exterior.coords)  # get points of a shape
-        stop = len(points) - 1  # define where to stop
-        for i in np.arange(
-            len(points)
-        ):  # for every point, calculate angle and add 1 if True angle
-            if i == 0:
-                continue
-            elif i == stop:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[1])
-                ang = _angle(a, b, c)
-
-                if ang <= 175:
-                    angles.append(ang)
-                elif _angle(a, b, c) >= 185:
-                    angles.append(ang)
-                else:
-                    continue
 
-            else:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[i + 1])
-                ang = _angle(a, b, c)
-
-                if _angle(a, b, c) <= 175:
-                    angles.append(ang)
-                elif _angle(a, b, c) >= 185:
-                    angles.append(ang)
-                else:
+    def __init__(self, gdf):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
+
+        def _angle(a, b, c):
+            ba = a - b
+            bc = c - b
+
+            cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
+            angle = np.degrees(np.arccos(cosine_angle))
+
+            return angle
+
+        # fill new column with the value of area, iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            angles = []
+            points = list(row["geometry"].exterior.coords)  # get points of a shape
+            stop = len(points) - 1  # define where to stop
+            for i in np.arange(
+                len(points)
+            ):  # for every point, calculate angle and add 1 if True angle
+                if i == 0:
                     continue
-        deviations = []
-        for i in angles:
-            dev = abs(90 - i)
-            deviations.append(dev)
-        results_list.append(np.mean(deviations))
-
-    series = pd.Series(results_list, index=gdf.index)
-
-    print("Squareness calculated.")
-    return series
-
+                elif i == stop:
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[1])
+                    ang = _angle(a, b, c)
+
+                    if ang <= 175:
+                        angles.append(ang)
+                    elif _angle(a, b, c) >= 185:
+                        angles.append(ang)
+                    else:
+                        continue
 
-def equivalent_rectangular_index(gdf, areas=None, perimeters=None):
+                else:
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[i + 1])
+                    ang = _angle(a, b, c)
+
+                    if _angle(a, b, c) <= 175:
+                        angles.append(ang)
+                    elif _angle(a, b, c) >= 185:
+                        angles.append(ang)
+                    else:
+                        continue
+            deviations = []
+            for i in angles:
+                dev = abs(90 - i)
+                deviations.append(dev)
+            results_list.append(np.mean(deviations))
+
+        self.s = pd.Series(results_list, index=gdf.index)
+
+
+class EquivalentRectangularIndex:
     """
     Calculates equivalent rectangular index of each object in given geoDataFrame.
 
@@ -941,10 +930,16 @@ def equivalent_rectangular_index(gdf, areas=None, perimeters=None):
         the name of the dataframe column, np.array, or pd.Series where is stored perimeter value. If set to None, function will calculate perimeters
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    eri : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
+    perimeters : Series
+        Series containing used perimeter values
 
     References
     ---------
@@ -956,46 +951,45 @@ def equivalent_rectangular_index(gdf, areas=None, perimeters=None):
 
     Examples
     --------
-    >>> buildings_df['eri'] = momepy.equivalent_rectangular_index(buildings_df, 'area', 'peri')
-    Calculating equivalent rectangular index...
+    >>> buildings_df['eri'] = momepy.EquivalentRectangularIndex(buildings_df, 'area', 'peri').eri
     100%|██████████| 144/144 [00:00<00:00, 895.57it/s]
-    Equivalent rectangular index calculated.
     >>> buildings_df['eri'][0]
     0.7879229963118455
     """
-    # define empty list for results
-    results_list = []
-    print("Calculating equivalent rectangular index...")
-    gdf = gdf.copy()
-
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+
+    def __init__(self, gdf, areas=None, perimeters=None):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
+        gdf = gdf.copy()
+
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
-    if areas is None:
-        gdf["mm_a"] = gdf.geometry.area
-        areas = "mm_a"
-    else:
-        if not isinstance(areas, str):
-            gdf["mm_a"] = areas
+        else:
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+        self.perimeters = gdf[perimeters]
+        if areas is None:
+            gdf["mm_a"] = gdf.geometry.area
             areas = "mm_a"
-    # fill new column with the value of area, iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        bbox = row["geometry"].minimum_rotated_rectangle
-        results_list.append(
-            math.sqrt(row[areas] / bbox.area) * (bbox.length / row[perimeters])
-        )
-
-    series = pd.Series(results_list, index=gdf.index)
+        else:
+            if not isinstance(areas, str):
+                gdf["mm_a"] = areas
+                areas = "mm_a"
+        self.areas = gdf[areas]
+        # fill new column with the value of area, iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            bbox = row["geometry"].minimum_rotated_rectangle
+            results_list.append(
+                math.sqrt(row[areas] / bbox.area) * (bbox.length / row[perimeters])
+            )
 
-    print("Equivalent rectangular index calculated.")
-    return series
+        self.eri = pd.Series(results_list, index=gdf.index)
 
 
-def elongation(gdf):
+class Elongation:
     """
     Calculates elongation of object seen as elongation of its minimum bounding rectangle.
 
@@ -1007,10 +1001,12 @@ def elongation(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    e : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     References
     ---------
@@ -1020,47 +1016,44 @@ def elongation(gdf):
 
     Examples
     --------
-    >>> buildings_df['elongation'] = momepy.elongation(buildings_df)
-    Calculating elongation...
+    >>> buildings_df['elongation'] = momepy.Elongation(buildings_df).e
     100%|██████████| 144/144 [00:00<00:00, 1032.62it/s]
-    Elongation calculated.
     >>> buildings_df['elongation'][0]
     0.9082437463675544
     """
-    # define empty list for results
-    results_list = []
-    print("Calculating elongation...")
-
-    # fill new column with the value of area, iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        bbox = row["geometry"].minimum_rotated_rectangle
-        a = bbox.area
-        p = bbox.length
-        cond1 = p ** 2
-        cond2 = 16 * a
-        if cond1 >= cond2:
-            sqrt = cond1 - cond2
-        else:
-            sqrt = 0
-
-        # calculate both width/length and length/width
-        elo1 = ((p - math.sqrt(sqrt)) / 4) / ((p / 2) - ((p - math.sqrt(sqrt)) / 4))
-        elo2 = ((p + math.sqrt(sqrt)) / 4) / ((p / 2) - ((p + math.sqrt(sqrt)) / 4))
-        # use the smaller one (e.g. shorter/longer)
-        if elo1 <= elo2:
-            elo = elo1
-        else:
-            elo = elo2
 
-        results_list.append(elo)
+    def __init__(self, gdf):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
+
+        # fill new column with the value of area, iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            bbox = row["geometry"].minimum_rotated_rectangle
+            a = bbox.area
+            p = bbox.length
+            cond1 = p ** 2
+            cond2 = 16 * a
+            if cond1 >= cond2:
+                sqrt = cond1 - cond2
+            else:
+                sqrt = 0
+
+            # calculate both width/length and length/width
+            elo1 = ((p - math.sqrt(sqrt)) / 4) / ((p / 2) - ((p - math.sqrt(sqrt)) / 4))
+            elo2 = ((p + math.sqrt(sqrt)) / 4) / ((p / 2) - ((p + math.sqrt(sqrt)) / 4))
+            # use the smaller one (e.g. shorter/longer)
+            if elo1 <= elo2:
+                elo = elo1
+            else:
+                elo = elo2
 
-    series = pd.Series(results_list, index=gdf.index)
+            results_list.append(elo)
 
-    print("Elongation calculated.")
-    return series
+        self.e = pd.Series(results_list, index=gdf.index)
 
 
-def centroid_corners(gdf):
+class CentroidCorners:
     """
     Calculates mean distance centroid - corners and st. deviation.
 
@@ -1072,14 +1065,14 @@ def centroid_corners(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects
 
-    Returns
-    -------
-    means, st_devs : tuple
-
-    means : Series
+    Attributes
+    ----------
+    mean : Series
         Series containing mean distance values.
-    st_devs : Series
+    std : Series
         Series containing standard deviation values.
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     References
     ----------
@@ -1089,90 +1082,87 @@ def centroid_corners(gdf):
 
     Examples
     --------
-    >>> means, deviations = momepy.centroid_corners(buildings_df)
-    Calculating distance centroid - corner...
+    >>> ccd = momepy.CentroidCorners(buildings_df)
     100%|██████████| 144/144 [00:00<00:00, 846.58it/s]
-    Distances centroid - corner calculated.
-    >>> buildings_df['ccd_means'] = means
-    >>> buildings_df['ccd_stdev'] = deviations
+    >>> buildings_df['ccd_means'] = ccd.means
+    >>> buildings_df['ccd_stdev'] = ccd.std
     >>> buildings_df['ccd_means'][0]
     15.961531913184833
     >>> buildings_df['ccd_stdev'][0]
     3.0810634305400177
     """
-    # define empty list for results
-    results_list = []
-    results_list_sd = []
-    print("Calculating distance centroid - corner...")
-
-    # calculate angle between points, return true or false if real corner
-    def true_angle(a, b, c):
-        ba = a - b
-        bc = c - b
-
-        cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
-        angle = np.arccos(cosine_angle)
-
-        if np.degrees(angle) <= 170:
-            return True
-        if np.degrees(angle) >= 190:
-            return True
-        return False
-
-    # iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        distances = []  # set empty list of distances
-        centroid = row["geometry"].centroid  # define centroid
-        points = list(row["geometry"].exterior.coords)  # get points of a shape
-        stop = len(points) - 1  # define where to stop
-        for i in np.arange(
-            len(points)
-        ):  # for every point, calculate angle and add 1 if True angle
-            if i == 0:
-                continue
-            elif i == stop:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[1])
-                p = Point(points[i])
-
-                if true_angle(a, b, c) is True:
-                    distance = centroid.distance(
-                        p
-                    )  # calculate distance point - centroid
-                    distances.append(distance)  # add distance to the list
-                else:
-                    continue
 
-            else:
-                a = np.asarray(points[i - 1])
-                b = np.asarray(points[i])
-                c = np.asarray(points[i + 1])
-                p = Point(points[i])
-
-                if true_angle(a, b, c) is True:
-                    distance = centroid.distance(p)
-                    distances.append(distance)
-                else:
+    def __init__(self, gdf):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
+        results_list_sd = []
+
+        # calculate angle between points, return true or false if real corner
+        def true_angle(a, b, c):
+            ba = a - b
+            bc = c - b
+
+            cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
+            angle = np.arccos(cosine_angle)
+
+            if np.degrees(angle) <= 170:
+                return True
+            if np.degrees(angle) >= 190:
+                return True
+            return False
+
+        # iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            distances = []  # set empty list of distances
+            centroid = row["geometry"].centroid  # define centroid
+            points = list(row["geometry"].exterior.coords)  # get points of a shape
+            stop = len(points) - 1  # define where to stop
+            for i in np.arange(
+                len(points)
+            ):  # for every point, calculate angle and add 1 if True angle
+                if i == 0:
                     continue
-        if not distances:  # circular buildings
-            from momepy.dimension import _longest_axis
-
-            results_list.append(
-                _longest_axis(row["geometry"].convex_hull.exterior.coords) / 2
-            )
-            results_list_sd.append(0)
-        else:
-            results_list.append(np.mean(distances))  # calculate mean
-            results_list_sd.append(np.std(distances))  # calculate st.dev
-    means = pd.Series(results_list, index=gdf.index)
-    st_devs = pd.Series(results_list_sd, index=gdf.index)
+                elif i == stop:
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[1])
+                    p = Point(points[i])
+
+                    if true_angle(a, b, c) is True:
+                        distance = centroid.distance(
+                            p
+                        )  # calculate distance point - centroid
+                        distances.append(distance)  # add distance to the list
+                    else:
+                        continue
 
-    print("Distances centroid - corner calculated.")
-    return means, st_devs
+                else:
+                    a = np.asarray(points[i - 1])
+                    b = np.asarray(points[i])
+                    c = np.asarray(points[i + 1])
+                    p = Point(points[i])
+
+                    if true_angle(a, b, c) is True:
+                        distance = centroid.distance(p)
+                        distances.append(distance)
+                    else:
+                        continue
+            if not distances:  # circular buildings
+                from momepy.dimension import _longest_axis
+
+                results_list.append(
+                    _longest_axis(row["geometry"].convex_hull.exterior.coords) / 2
+                )
+                results_list_sd.append(0)
+            else:
+                results_list.append(np.mean(distances))  # calculate mean
+                results_list_sd.append(np.std(distances))  # calculate st.dev
+        self.mean = pd.Series(results_list, index=gdf.index)
+        self.std = pd.Series(results_list_sd, index=gdf.index)
 
 
-def linearity(gdf):
+class Linearity:
     """
     Calculates linearity of each LineString object in given geoDataFrame.
 
@@ -1184,10 +1174,12 @@ def linearity(gdf):
     gdf : GeoDataFrame
         GeoDataFrame containing objects
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    linearity : Series
+        Series containing mean distance values.
+    gdf : GeoDataFrame
+        original GeoDataFrame
 
     References
     ---------
@@ -1197,31 +1189,28 @@ def linearity(gdf):
 
     Examples
     --------
-    >>> streets_df['linearity'] = momepy.linearity(streets_df)
-    Calculating linearity...
+    >>> streets_df['linearity'] = momepy.Linearity(streets_df).linearity
     100%|██████████| 33/33 [00:00<00:00, 1737.64it/s]
-    Linearity calculated.
     >>> streets_df['linearity'][0]
     1.0
     """
-    # define empty list for results
-    results_list = []
-    print("Calculating linearity...")
-
-    # fill new column with the value of area, iterating over rows one by one
-    for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
-        euclidean = Point(row["geometry"].coords[0]).distance(
-            Point(row["geometry"].coords[-1])
-        )
-        results_list.append(euclidean / row["geometry"].length)
 
-    series = pd.Series(results_list, index=gdf.index)
+    def __init__(self, gdf):
+        self.gdf = gdf
+        # define empty list for results
+        results_list = []
 
-    print("Linearity calculated.")
-    return series
+        # fill new column with the value of area, iterating over rows one by one
+        for index, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
+            euclidean = Point(row["geometry"].coords[0]).distance(
+                Point(row["geometry"].coords[-1])
+            )
+            results_list.append(euclidean / row["geometry"].length)
+
+        self.linearity = pd.Series(results_list, index=gdf.index)
 
 
-def compactness_weighted_axis(gdf, areas=None, perimeters=None, longest_axis=None):
+class CompactnessWeightedAxis:
     """
     Calculates compactness-weighted axis of each object in given geoDataFrame.
 
@@ -1244,54 +1233,55 @@ def compactness_weighted_axis(gdf, areas=None, perimeters=None, longest_axis=Non
         the name of the dataframe column, np.array, or pd.Series where is stored longest axis length value. If set to None, function will calculate it
         during the process without saving them separately.
 
-    Returns
-    -------
-    Series
-        Series containing resulting values.
+    Attributes
+    ----------
+    cwa : Series
+        Series containing resulting values
+    gdf : GeoDataFrame
+        original GeoDataFrame
+    areas : Series
+        Series containing used area values
+    longest_axis : Series
+        Series containing used area values
+    perimeters : Series
+        Series containing used area values
+
 
     Examples
     --------
-    >>> blocks_df['cwa'] = mm.compactness_weighted_axis(blocks_df)
-    Calculating compactness-weighted axis...
-    Compactness-weighted axis calculated.
+    >>> blocks_df['cwa'] = mm.CompactnessWeightedAxis(blocks_df).cwa
     """
 
-    print("Calculating compactness-weighted axis...")
-    gdf = gdf.copy()
+    def __init__(self, gdf, areas=None, perimeters=None, longest_axis=None):
+        self.gdf = gdf
+        gdf = gdf.copy()
 
-    if perimeters is None:
-        gdf["mm_p"] = gdf.geometry.length
-        perimeters = "mm_p"
-    else:
-        if not isinstance(perimeters, str):
-            gdf["mm_p"] = perimeters
+        if perimeters is None:
+            gdf["mm_p"] = gdf.geometry.length
             perimeters = "mm_p"
-
-    if longest_axis is None:
-        from .dimension import longest_axis_length
-
-        gdf["mm_la"] = longest_axis_length(gdf)
-        longest_axis = "mm_la"
-    else:
-        if not isinstance(longest_axis, str):
-            gdf["mm_la"] = longest_axis
+        else:
+            if not isinstance(perimeters, str):
+                gdf["mm_p"] = perimeters
+                perimeters = "mm_p"
+        self.perimeters = gdf[perimeters]
+        if longest_axis is None:
+            from .dimension import LongestAxisLength
+
+            gdf["mm_la"] = LongestAxisLength(gdf).lal
             longest_axis = "mm_la"
-
-    if areas is None:
-        series = gdf.apply(
-            lambda row: row[longest_axis]
-            * ((4 / math.pi) - (16 * row.geometry.area) / ((row[perimeters]) ** 2)),
-            axis=1,
-        )
-    else:
+        else:
+            if not isinstance(longest_axis, str):
+                gdf["mm_la"] = longest_axis
+                longest_axis = "mm_la"
+        self.longest_axis = gdf[longest_axis]
+        if areas is None:
+            areas = gdf.geometry.area
         if not isinstance(areas, str):
             gdf["mm_a"] = areas
             areas = "mm_a"
-        series = gdf.apply(
+        self.areas = gdf[areas]
+        self.cwa = gdf.apply(
             lambda row: row[longest_axis]
             * ((4 / math.pi) - (16 * row[areas]) / ((row[perimeters]) ** 2)),
             axis=1,
         )
-
-    print("Compactness-weighted axis calculated.")
-    return series
diff --git a/momepy/utils.py b/momepy/utils.py
index c7eb871..b2cc334 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -11,7 +11,7 @@ import operator
 import shapely
 import math
 
-from .shape import circular_compactness
+from .shape import CircularCompactness
 
 
 def unique_id(objects):
@@ -343,10 +343,8 @@ def limit_range(vals, rng):
         rng = sorted(rng)
         lower = np.percentile(vals, rng[0], interpolation="nearest")
         higher = np.percentile(vals, rng[1], interpolation="nearest")
-        for x in vals:
-            if x >= lower and x <= higher:
-                limited.append(x)
-        return limited
+        limited = [x for x in vals if x >= lower and x <= higher]
+        return np.array(limited)
     return vals
 
 
@@ -399,7 +397,7 @@ def preprocess(buildings, size=30, compactness=True, islands=True):
         blg["neighbors"] = sw.neighbors
         blg["neighbors"] = blg["neighbors"].map(sw.neighbors)
         blg["n_count"] = blg.apply(lambda row: len(row.neighbors), axis=1)
-        blg["circu"] = circular_compactness(blg)
+        blg["circu"] = CircularCompactness(blg).cc
 
         # idetify those smaller than x with only one neighbor and attaches it to it.
         join = {}


Test patch:
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index fc2435d..51e0575 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -18,139 +18,128 @@ class TestDimensions:
         self.df_tessellation = gpd.read_file(test_file_path, layer="tessellation")
         self.df_buildings["height"] = np.linspace(10.0, 30.0, 144)
 
-    def test_area(self):
-        self.df_buildings["area"] = mm.area(self.df_buildings)
+    def test_Area(self):
+        self.df_buildings["area"] = mm.Area(self.df_buildings).area
         check = self.df_buildings.geometry[0].area
         assert self.df_buildings["area"][0] == check
 
-    def test_perimeter(self):
-        self.df_buildings["perimeter"] = mm.perimeter(self.df_buildings)
+    def test_Perimeter(self):
+        self.df_buildings["perimeter"] = mm.Perimeter(self.df_buildings).perimeter
         check = self.df_buildings.geometry[0].length
         assert self.df_buildings["perimeter"][0] == check
 
-    def test_volume(self):
+    def test_Volume(self):
         self.df_buildings["area"] = self.df_buildings.geometry.area
-        self.df_buildings["volume"] = mm.volume(self.df_buildings, "height", "area")
+        self.df_buildings["volume"] = mm.Volume(
+            self.df_buildings, "height", "area"
+        ).volume
         check = self.df_buildings.geometry[0].area * self.df_buildings.height[0]
         assert self.df_buildings["volume"][0] == check
 
-    def test_volume_aray(self):
         area = self.df_buildings.geometry.area
         height = np.linspace(10.0, 30.0, 144)
-        self.df_buildings["volume"] = mm.volume(self.df_buildings, height, area)
+        self.df_buildings["volume"] = mm.Volume(self.df_buildings, height, area).volume
         check = self.df_buildings.geometry[0].area * self.df_buildings.height[0]
         assert self.df_buildings["volume"][0] == check
 
-    def test_volume_no_area(self):
-        self.df_buildings["volume"] = mm.volume(self.df_buildings, "height")
+        self.df_buildings["volume"] = mm.Volume(self.df_buildings, "height").volume
         check = self.df_buildings.geometry[0].area * self.df_buildings.height[0]
         assert self.df_buildings["volume"][0] == check
 
-    def test_volume_missing_col(self):
         with pytest.raises(KeyError):
-            self.df_buildings["volume"] = mm.volume(self.df_buildings, "height", "area")
+            self.df_buildings["volume"] = mm.Volume(
+                self.df_buildings, "height", "nonexistent"
+            )
 
-    def test_floor_area(self):
+    def test_FloorArea(self):
         self.df_buildings["area"] = self.df_buildings.geometry.area
-        self.df_buildings["floor_area"] = mm.floor_area(
+        self.df_buildings["floor_area"] = mm.FloorArea(
             self.df_buildings, "height", "area"
-        )
+        ).fa
         check = self.df_buildings.geometry[0].area * (self.df_buildings.height[0] // 3)
         assert self.df_buildings["floor_area"][0] == check
 
-    def test_floor_area_array(self):
         area = self.df_buildings.geometry.area
         height = np.linspace(10.0, 30.0, 144)
-        self.df_buildings["floor_area"] = mm.floor_area(self.df_buildings, height, area)
-        check = self.df_buildings.geometry[0].area * (self.df_buildings.height[0] // 3)
+        self.df_buildings["floor_area"] = mm.FloorArea(
+            self.df_buildings, height, area
+        ).fa
         assert self.df_buildings["floor_area"][0] == check
 
-    def test_floor_area_no_area(self):
-        self.df_buildings["floor_area"] = mm.floor_area(self.df_buildings, "height")
-        check = self.df_buildings.geometry[0].area * (self.df_buildings.height[0] // 3)
+        self.df_buildings["floor_area"] = mm.FloorArea(self.df_buildings, "height").fa
         assert self.df_buildings["floor_area"][0] == check
 
-    def test_floor_area_missing_col(self):
         with pytest.raises(KeyError):
-            self.df_buildings["floor_area"] = mm.floor_area(
-                self.df_buildings, "height", "area"
+            self.df_buildings["floor_area"] = mm.FloorArea(
+                self.df_buildings, "height", "nonexistent"
             )
 
-    def test_courtyard_area(self):
+    def test_CourtyardArea(self):
         self.df_buildings["area"] = self.df_buildings.geometry.area
-        self.df_buildings["courtyard_area"] = mm.courtyard_area(
+        self.df_buildings["courtyard_area"] = mm.CourtyardArea(
             self.df_buildings, "area"
-        )
+        ).ca
         check = (
             Polygon(self.df_buildings.geometry[80].exterior).area
             - self.df_buildings.geometry[80].area
         )
         assert self.df_buildings["courtyard_area"][80] == check
 
-    def test_courtyard_area_array(self):
         area = self.df_buildings.geometry.area
-        self.df_buildings["courtyard_area"] = mm.courtyard_area(self.df_buildings, area)
-        check = (
-            Polygon(self.df_buildings.geometry[80].exterior).area
-            - self.df_buildings.geometry[80].area
-        )
+        self.df_buildings["courtyard_area"] = mm.CourtyardArea(
+            self.df_buildings, area
+        ).ca
         assert self.df_buildings["courtyard_area"][80] == check
 
-    def test_courtyard_no_area(self):
-        self.df_buildings["courtyard_area"] = mm.courtyard_area(self.df_buildings)
-        check = (
-            Polygon(self.df_buildings.geometry[80].exterior).area
-            - self.df_buildings.geometry[80].area
-        )
+        self.df_buildings["courtyard_area"] = mm.CourtyardArea(self.df_buildings).ca
         assert self.df_buildings["courtyard_area"][80] == check
 
-    def test_courtyard_missing_col(self):
         with pytest.raises(KeyError):
-            self.df_buildings["courtyard_area"] = mm.floor_area(
-                self.df_buildings, "area"
+            self.df_buildings["courtyard_area"] = mm.CourtyardArea(
+                self.df_buildings, "nonexistent"
             )
 
-    def test_longest_axis_length(self):
-        self.df_buildings["long_axis"] = mm.longest_axis_length(self.df_buildings)
+    def test_LongestAxisLength(self):
+        self.df_buildings["long_axis"] = mm.LongestAxisLength(self.df_buildings).lal
         check = (
             _make_circle(self.df_buildings.geometry[0].convex_hull.exterior.coords)[2]
             * 2
         )
         assert self.df_buildings["long_axis"][0] == check
 
-    def test_average_character(self):
+    def test_AverageCharacter(self):
         spatial_weights = sw_high(k=3, gdf=self.df_tessellation, ids="uID")
         self.df_tessellation["area"] = area = self.df_tessellation.geometry.area
-        self.df_tessellation["mesh_ar"] = mm.average_character(
+        self.df_tessellation["mesh_ar"] = mm.AverageCharacter(
             self.df_tessellation,
             values="area",
             spatial_weights=spatial_weights,
             unique_id="uID",
             mode="mode",
-        )
-        self.df_tessellation["mesh_array"] = mm.average_character(
+        ).ac
+        self.df_tessellation["mesh_array"] = mm.AverageCharacter(
             self.df_tessellation,
             values=area,
             spatial_weights=spatial_weights,
             unique_id="uID",
             mode="median",
-        )
-        self.df_tessellation["mesh_id"] = mm.average_character(
+        ).ac
+        self.df_tessellation["mesh_id"] = mm.AverageCharacter(
             self.df_tessellation,
             spatial_weights=spatial_weights,
             values="area",
             rng=(10, 90),
             unique_id="uID",
-        )
-        self.df_tessellation["mesh_iq"] = mm.average_character(
+        ).ac
+        self.df_tessellation["mesh_iq"] = mm.AverageCharacter(
             self.df_tessellation,
             spatial_weights=spatial_weights,
             values="area",
             rng=(25, 75),
             unique_id="uID",
-        )
+        ).ac
         with pytest.raises(ValueError):
-            self.df_tessellation["mesh_ar"] = mm.average_character(
+            self.df_tessellation["mesh_ar"] = mm.AverageCharacter(
                 self.df_tessellation,
                 values="area",
                 spatial_weights=spatial_weights,
@@ -162,62 +151,59 @@ class TestDimensions:
         assert self.df_tessellation["mesh_id"][38] == approx(2250.224, rel=1e-3)
         assert self.df_tessellation["mesh_iq"][38] == approx(2118.609, rel=1e-3)
 
-    def test_street_profile(self):
-        results = mm.street_profile(
-            self.df_streets, self.df_buildings, heights="height"
-        )
-        assert results["widths"][0] == 46.7290758769204
-        assert results["width_deviations"][0] == 0.5690458650581023
-        assert results["heights"][0] == 19.454545454545457
-        assert results["profile"][0] == 0.4163263469148574
-        assert results["openness"][0] == 0.9107142857142857
-        assert results["heights_deviations"][0] == 7.720786443287433
-
-    def test_street_profile_array(self):
+    def test_StreetProfile(self):
+        results = mm.StreetProfile(self.df_streets, self.df_buildings, heights="height")
+        assert results.w[0] == 46.7290758769204
+        assert results.wd[0] == 0.5690458650581023
+        assert results.h[0] == 19.454545454545457
+        assert results.p[0] == 0.4163263469148574
+        assert results.o[0] == 0.9107142857142857
+        assert results.hd[0] == 7.720786443287433
+
         height = np.linspace(10.0, 30.0, 144)
-        results = mm.street_profile(
+        results2 = mm.StreetProfile(
             self.df_streets, self.df_buildings, heights=height, tick_length=100
         )
-        assert results["widths"][0] == 67.563796073073
-        assert results["width_deviations"][0] == 8.791875291865827
-        assert results["heights"][0] == 23.756643356643362
-        assert results["profile"][0] == 0.3516179483306353
-        assert results["openness"][0] == 0.5535714285714286
-        assert results["heights_deviations"][0] == 5.526848034418866
-
-    def test_weighted_character_sw(self):
+        assert results2.w[0] == 67.563796073073
+        assert results2.wd[0] == 8.791875291865827
+        assert results2.h[0] == 23.756643356643362
+        assert results2.p[0] == 0.3516179483306353
+        assert results2.o[0] == 0.5535714285714286
+        assert results2.hd[0] == 5.526848034418866
+
+    def test_WeightedCharacter(self):
         sw = sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        weighted = mm.weighted_character(self.df_buildings, "height", sw, "uID")
+        weighted = mm.WeightedCharacter(self.df_buildings, "height", sw, "uID").wc
         assert weighted[38] == approx(18.301, rel=1e-3)
 
-    def test_weighted_character_area(self):
         self.df_buildings["area"] = self.df_buildings.geometry.area
         sw = sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        weighted = mm.weighted_character(self.df_buildings, "height", sw, "uID", "area")
+        weighted = mm.WeightedCharacter(
+            self.df_buildings, "height", sw, "uID", "area"
+        ).wc
         assert weighted[38] == approx(18.301, rel=1e-3)
 
-    def test_weighted_character_array(self):
         area = self.df_buildings.geometry.area
         sw = sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        weighted = mm.weighted_character(
+        weighted = mm.WeightedCharacter(
             self.df_buildings, self.df_buildings.height, sw, "uID", area
-        )
+        ).wc
         assert weighted[38] == approx(18.301, rel=1e-3)
 
-    def test_covered_area(self):
+    def test_CoveredArea(self):
         sw = sw_high(gdf=self.df_tessellation, k=1, ids="uID")
-        covered_sw = mm.covered_area(self.df_tessellation, sw, "uID")
+        covered_sw = mm.CoveredArea(self.df_tessellation, sw, "uID").ca
         assert covered_sw[0] == approx(24115.667, rel=1e-3)
 
-    def test_wall(self):
+    def test_PerimeterWall(self):
         sw = sw_high(gdf=self.df_buildings, k=1)
-        wall = mm.wall(self.df_buildings)
-        wall_sw = mm.wall(self.df_buildings, sw)
+        wall = mm.PerimeterWall(self.df_buildings).wall
+        wall_sw = mm.PerimeterWall(self.df_buildings, sw).wall
         assert wall[0] == wall_sw[0]
         assert wall[0] == 137.2106961418436
 
-    def test_segments_length(self):
-        absol = mm.segments_length(self.df_streets)
-        mean = mm.segments_length(self.df_streets, mean=True)
+    def test_SegmentsLength(self):
+        absol = mm.SegmentsLength(self.df_streets).sl
+        mean = mm.SegmentsLength(self.df_streets, mean=True).sl
         assert max(absol) == 1907.502238338006
         assert max(mean) == 249.5698434867373
diff --git a/tests/test_distribution.py b/tests/test_distribution.py
index bb5cf86..25a7449 100644
--- a/tests/test_distribution.py
+++ b/tests/test_distribution.py
@@ -14,60 +14,60 @@ class TestDistribution:
         self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.df_tessellation = gpd.read_file(test_file_path, layer="tessellation")
         self.df_buildings["height"] = np.linspace(10.0, 30.0, 144)
-        self.df_buildings["volume"] = mm.volume(self.df_buildings, "height")
+        self.df_buildings["volume"] = mm.Volume(self.df_buildings, "height").volume
         self.df_streets["nID"] = mm.unique_id(self.df_streets)
         self.df_buildings["nID"] = mm.get_network_id(
             self.df_buildings, self.df_streets, "uID", "nID"
         )
 
-    def test_orientation(self):
-        self.df_buildings["orient"] = mm.orientation(self.df_buildings)
+    def test_Orientation(self):
+        self.df_buildings["orient"] = mm.Orientation(self.df_buildings).o
         check = 41.05146788287027
         assert self.df_buildings["orient"][0] == check
 
-    def test_shared_walls_ratio(self):
-        self.df_buildings["swr"] = mm.shared_walls_ratio(self.df_buildings, "uID")
-        self.df_buildings["swr_uid"] = mm.shared_walls_ratio(
+    def test_SharedWallsRatio(self):
+        self.df_buildings["swr"] = mm.SharedWallsRatio(self.df_buildings, "uID").swr
+        self.df_buildings["swr_uid"] = mm.SharedWallsRatio(
             self.df_buildings, range(len(self.df_buildings))
-        )
-        self.df_buildings["swr_array"] = mm.shared_walls_ratio(
+        ).swr
+        self.df_buildings["swr_array"] = mm.SharedWallsRatio(
             self.df_buildings, "uID", self.df_buildings.geometry.length
-        )
+        ).swr
         check = 0.3424804411228673
         assert self.df_buildings["swr"][10] == check
         assert self.df_buildings["swr_uid"][10] == check
         assert self.df_buildings["swr_array"][10] == check
 
-    def test_street_alignment(self):
-        self.df_buildings["orient"] = orient = mm.orientation(self.df_buildings)
-        self.df_buildings["street_alignment"] = mm.street_alignment(
+    def test_StreetAlignment(self):
+        self.df_buildings["orient"] = orient = mm.Orientation(self.df_buildings).o
+        self.df_buildings["street_alignment"] = mm.StreetAlignment(
             self.df_buildings, self.df_streets, "orient", network_id="nID"
-        )
-        self.df_buildings["street_alignment2"] = mm.street_alignment(
+        ).sa
+        self.df_buildings["street_alignment2"] = mm.StreetAlignment(
             self.df_buildings,
             self.df_streets,
             "orient",
             left_network_id="nID",
             right_network_id="nID",
-        )
-        self.df_buildings["street_a_arr"] = mm.street_alignment(
+        ).sa
+        self.df_buildings["street_a_arr"] = mm.StreetAlignment(
             self.df_buildings,
             self.df_streets,
             orient,
             left_network_id=self.df_buildings["nID"],
             right_network_id=self.df_streets["nID"],
-        )
+        ).sa
 
         with pytest.raises(ValueError):
-            self.df_buildings["street_alignment"] = mm.street_alignment(
+            self.df_buildings["street_alignment"] = mm.StreetAlignment(
                 self.df_buildings, self.df_streets, "orient"
             )
         with pytest.raises(ValueError):
-            self.df_buildings["street_alignment"] = mm.street_alignment(
+            self.df_buildings["street_alignment"] = mm.StreetAlignment(
                 self.df_buildings, self.df_streets, "orient", left_network_id="nID"
             )
         with pytest.raises(ValueError):
-            self.df_buildings["street_alignment"] = mm.street_alignment(
+            self.df_buildings["street_alignment"] = mm.StreetAlignment(
                 self.df_buildings, self.df_streets, "orient", right_network_id="nID"
             )
         check = 0.29073888476702336
@@ -75,15 +75,17 @@ class TestDistribution:
         assert self.df_buildings["street_alignment2"][0] == check
         assert self.df_buildings["street_a_arr"][0] == check
 
-    def test_cell_alignment(self):
-        self.df_buildings["orient"] = blgori = mm.orientation(self.df_buildings)
-        self.df_tessellation["orient"] = tessori = mm.orientation(self.df_tessellation)
-        self.df_buildings["c_align"] = mm.cell_alignment(
+    def test_CellAlignment(self):
+        self.df_buildings["orient"] = blgori = mm.Orientation(self.df_buildings).o
+        self.df_tessellation["orient"] = tessori = mm.Orientation(
+            self.df_tessellation
+        ).o
+        self.df_buildings["c_align"] = mm.CellAlignment(
             self.df_buildings, self.df_tessellation, "orient", "orient", "uID", "uID"
-        )
-        self.df_buildings["c_align_array"] = mm.cell_alignment(
+        ).ca
+        self.df_buildings["c_align_array"] = mm.CellAlignment(
             self.df_buildings, self.df_tessellation, blgori, tessori, "uID", "uID"
-        )
+        ).ca
         check = abs(
             self.df_buildings["orient"][0]
             - self.df_tessellation[
@@ -92,64 +94,64 @@ class TestDistribution:
         )
         assert self.df_buildings["c_align"][0] == check
 
-    def test_alignment(self):
-        self.df_buildings["orient"] = mm.orientation(self.df_buildings)
+    def test_Alignment(self):
+        self.df_buildings["orient"] = mm.Orientation(self.df_buildings).o
         sw = Queen.from_dataframe(self.df_tessellation, ids="uID")
-        self.df_buildings["align_sw"] = mm.alignment(
+        self.df_buildings["align_sw"] = mm.Alignment(
             self.df_buildings, sw, "uID", self.df_buildings["orient"]
-        )
+        ).a
         assert self.df_buildings["align_sw"][0] == 18.299481296455237
 
-    def test_neighbour_distance(self):
+    def test_NeighborDistance(self):
         sw = Queen.from_dataframe(self.df_tessellation, ids="uID")
-        self.df_buildings["dist_sw"] = mm.neighbour_distance(
+        self.df_buildings["dist_sw"] = mm.NeighborDistance(
             self.df_buildings, sw, "uID"
-        )
+        ).nd
         check = 29.18589019096464
         assert self.df_buildings["dist_sw"][0] == check
 
-    def test_mean_interbuilding_distance(self):
+    def test_MeanInterbuildingDistance(self):
         sw = Queen.from_dataframe(self.df_tessellation, ids="uID")
         swh = mm.sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        self.df_buildings["m_dist_sw"] = mm.mean_interbuilding_distance(
+        self.df_buildings["m_dist_sw"] = mm.MeanInterbuildingDistance(
             self.df_buildings, sw, "uID", swh
-        )
-        self.df_buildings["m_dist"] = mm.mean_interbuilding_distance(
+        ).mid
+        self.df_buildings["m_dist"] = mm.MeanInterbuildingDistance(
             self.df_buildings, sw, "uID", order=3
-        )
+        ).mid
         check = 29.305457092042744
         assert self.df_buildings["m_dist_sw"][0] == check
         assert self.df_buildings["m_dist"][0] == check
 
-    def test_neighbouring_street_orientation_deviation(self):
-        self.df_streets["dev"] = mm.neighbouring_street_orientation_deviation(
+    def test_NeighboringStreetOrientationDeviation(self):
+        self.df_streets["dev"] = mm.NeighboringStreetOrientationDeviation(
             self.df_streets
-        )
+        ).nsod
         check = 5.986848512501008
         assert self.df_streets["dev"].mean() == check
 
-    def test_building_adjacency(self):
+    def test_BuildingAdjacencyy(self):
         sw = Queen.from_dataframe(self.df_buildings, ids="uID")
         swh = mm.sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        self.df_buildings["adj_sw"] = mm.building_adjacency(
+        self.df_buildings["adj_sw"] = mm.BuildingAdjacency(
             self.df_buildings,
             spatial_weights=sw,
             unique_id="uID",
             spatial_weights_higher=swh,
-        )
-        self.df_buildings["adj_sw_none"] = mm.building_adjacency(
+        ).ba
+        self.df_buildings["adj_sw_none"] = mm.BuildingAdjacency(
             self.df_buildings, unique_id="uID", spatial_weights_higher=swh
-        )
+        ).ba
         check = 0.2613824113909074
         assert self.df_buildings["adj_sw"].mean() == check
         assert self.df_buildings["adj_sw_none"].mean() == check
 
-    def test_neighbours(self):
+    def test_Neighbors(self):
         sw = Queen.from_dataframe(self.df_tessellation, ids="uID")
-        self.df_tessellation["nei_sw"] = mm.neighbours(self.df_tessellation, sw, "uID")
-        self.df_tessellation["nei_wei"] = mm.neighbours(
+        self.df_tessellation["nei_sw"] = mm.Neighbors(self.df_tessellation, sw, "uID").n
+        self.df_tessellation["nei_wei"] = mm.Neighbors(
             self.df_tessellation, sw, "uID", weighted=True
-        )
+        ).n
         check = 5.180555555555555
         check_w = 0.029066398893536072
         assert self.df_tessellation["nei_sw"].mean() == check
diff --git a/tests/test_diversity.py b/tests/test_diversity.py
index 6480940..d432042 100644
--- a/tests/test_diversity.py
+++ b/tests/test_diversity.py
@@ -16,56 +16,58 @@ class TestDiversity:
         self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.df_tessellation = gpd.read_file(test_file_path, layer="tessellation")
         self.df_buildings["height"] = np.linspace(10.0, 30.0, 144)
-        self.df_tessellation["area"] = mm.area(self.df_tessellation)
+        self.df_tessellation["area"] = mm.Area(self.df_tessellation).area
         self.sw = sw_high(k=3, gdf=self.df_tessellation, ids="uID")
 
-    def test_rng(self):
-        full_sw = mm.rng(self.df_tessellation, "area", self.sw, "uID")
+    def test_Range(self):
+        full_sw = mm.Range(self.df_tessellation, "area", self.sw, "uID").r
         assert full_sw[0] == approx(8255.372, rel=1e-3)
         area = self.df_tessellation["area"]
-        full2 = mm.rng(self.df_tessellation, area, self.sw, "uID")
+        full2 = mm.Range(self.df_tessellation, area, self.sw, "uID").r
         assert full2[0] == approx(8255.372, rel=1e-3)
-        limit = mm.rng(self.df_tessellation, "area", self.sw, "uID", rng=(10, 90))
+        limit = mm.Range(self.df_tessellation, "area", self.sw, "uID", rng=(10, 90)).r
         assert limit[0] == approx(4122.139, rel=1e-3)
 
-    def test_theil(self):
-        full_sw = mm.theil(self.df_tessellation, "area", self.sw, "uID")
+    def test_Theil(self):
+        full_sw = mm.Theil(self.df_tessellation, "area", self.sw, "uID").t
         assert full_sw[0] == approx(0.25744684)
-        limit = mm.theil(
+        limit = mm.Theil(
             self.df_tessellation,
             self.df_tessellation.area,
             self.sw,
             "uID",
             rng=(10, 90),
-        )
+        ).t
         assert limit[0] == approx(0.1330295)
-        zeros = mm.theil(self.df_tessellation, np.zeros(len(self.df_tessellation)), self.sw, "uID")
+        zeros = mm.Theil(
+            self.df_tessellation, np.zeros(len(self.df_tessellation)), self.sw, "uID"
+        ).t
         assert zeros[0] == 0
 
-    def test_simpson(self):
-        ht_sw = mm.simpson(self.df_tessellation, "area", self.sw, "uID")
+    def test_Simpson(self):
+        ht_sw = mm.Simpson(self.df_tessellation, "area", self.sw, "uID").s
         assert ht_sw[0] == 0.385
-        quan_sw = mm.simpson(
+        quan_sw = mm.Simpson(
             self.df_tessellation,
             self.df_tessellation.area,
             self.sw,
             "uID",
             binning="quantiles",
             k=3,
-        )
+        ).s
         assert quan_sw[0] == 0.395
         with pytest.raises(ValueError):
-            ht_sw = mm.simpson(
+            ht_sw = mm.Simpson(
                 self.df_tessellation, "area", self.sw, "uID", binning="nonexistent"
             )
 
-    def test_gini(self):
-        full_sw = mm.gini(self.df_tessellation, "area", self.sw, "uID")
+    def test_Gini(self):
+        full_sw = mm.Gini(self.df_tessellation, "area", self.sw, "uID").g
         assert full_sw[0] == approx(0.3945388)
-        limit = mm.gini(self.df_tessellation, "area", self.sw, "uID", rng=(10, 90))
+        limit = mm.Gini(self.df_tessellation, "area", self.sw, "uID", rng=(10, 90)).g
         assert limit[0] == approx(0.28532814)
         self.df_tessellation["negative"] = (
             self.df_tessellation.area - self.df_tessellation.area.mean()
         )
         with pytest.raises(ValueError):
-            mm.gini(self.df_tessellation, "negative", self.sw, "uID")
+            mm.Gini(self.df_tessellation, "negative", self.sw, "uID").g
diff --git a/tests/test_intensity.py b/tests/test_intensity.py
index 19452eb..16eb9f4 100644
--- a/tests/test_intensity.py
+++ b/tests/test_intensity.py
@@ -19,7 +19,7 @@ class TestIntensity:
         self.df_buildings["height"] = np.linspace(10.0, 30.0, 144)
         self.df_tessellation["area"] = self.df_tessellation.geometry.area
         self.df_buildings["area"] = self.df_buildings.geometry.area
-        self.df_buildings["fl_area"] = mm.floor_area(self.df_buildings, "height")
+        self.df_buildings["fl_area"] = mm.FloorArea(self.df_buildings, "height").fa
         self.df_buildings["nID"] = mm.get_network_id(
             self.df_buildings, self.df_streets, "uID", "nID"
         )
@@ -27,36 +27,34 @@ class TestIntensity:
             self.df_tessellation, self.df_streets, self.df_buildings, "bID", "uID"
         )
 
-    def test_object_area_ratio(self):
-        car = mm.object_area_ratio(
+    def test_AreaRatio(self):
+        car = mm.AreaRatio(
             self.df_tessellation, self.df_buildings, "area", "area", "uID"
-        )
-        carlr = mm.object_area_ratio(
+        ).ar
+        carlr = mm.AreaRatio(
             self.df_tessellation,
             self.df_buildings,
             "area",
             "area",
             left_unique_id="uID",
             right_unique_id="uID",
-        )
+        ).ar
         check = 0.3206556897709747
         assert car.mean() == check
         assert carlr.mean() == check
-        far = mm.object_area_ratio(
+        far = mm.AreaRatio(
             self.df_tessellation,
             self.df_buildings,
             self.df_tessellation.area,
             self.df_buildings.fl_area,
             "uID",
-        )
+        ).ar
         check = 1.910949846262234
         assert far.mean() == check
         with pytest.raises(ValueError):
-            car = mm.object_area_ratio(
-                self.df_tessellation, self.df_buildings, "area", "area"
-            )
+            car = mm.AreaRatio(self.df_tessellation, self.df_buildings, "area", "area")
         with pytest.raises(ValueError):
-            car = mm.object_area_ratio(
+            car = mm.AreaRatio(
                 self.df_tessellation,
                 self.df_buildings,
                 "area",
@@ -64,7 +62,7 @@ class TestIntensity:
                 left_unique_id="uID",
             )
         with pytest.raises(ValueError):
-            car = mm.object_area_ratio(
+            car = mm.AreaRatio(
                 self.df_tessellation,
                 self.df_buildings,
                 "area",
@@ -72,84 +70,84 @@ class TestIntensity:
                 right_unique_id="uID",
             )
 
-    def test_elements_count(self):
-        eib = mm.elements_count(self.blocks, self.df_buildings, "bID", "bID")
-        weib = mm.elements_count(
-            self.blocks, self.df_buildings, "bID", "bID", weighted=True
-        )
-        weis = mm.elements_count(
+    def test_Count(self):
+        eib = mm.Count(self.blocks, self.df_buildings, "bID", "bID").c
+        weib = mm.Count(self.blocks, self.df_buildings, "bID", "bID", weighted=True).c
+        weis = mm.Count(
             self.df_streets, self.df_buildings, "nID", "nID", weighted=True
-        )
+        ).c
         check_eib = [13, 14, 8, 26, 24, 17, 23, 19]
         check_weib = 0.00040170607189453996
         assert eib.tolist() == check_eib
         assert weib.mean() == check_weib
         assert weis.mean() == 0.020524232642849215
 
-    def test_courtyards(self):
-        courtyards = mm.courtyards(self.df_buildings, "bID")
+    def test_Courtyards(self):
+        courtyards = mm.Courtyards(self.df_buildings, "bID").c
         sw = Queen.from_dataframe(self.df_buildings)
-        courtyards_wm = mm.courtyards(self.df_buildings, self.df_buildings.bID, sw)
+        courtyards_wm = mm.Courtyards(self.df_buildings, self.df_buildings.bID, sw).c
         check = 0.6805555555555556
         assert courtyards.mean() == check
         assert courtyards_wm.mean() == check
 
-    def test_blocks_count(self):
+    def test_BlocksCount(self):
         sw = mm.sw_high(k=5, gdf=self.df_tessellation, ids="uID")
-        count = mm.blocks_count(self.df_tessellation, "bID", sw, "uID")
-        count2 = mm.blocks_count(
+        count = mm.BlocksCount(self.df_tessellation, "bID", sw, "uID").bc
+        count2 = mm.BlocksCount(
             self.df_tessellation, self.df_tessellation.bID, sw, "uID"
-        )
-        unweigthed = mm.blocks_count(
+        ).bc
+        unweigthed = mm.BlocksCount(
             self.df_tessellation, "bID", sw, "uID", weighted=False
-        )
+        ).bc
         check = 3.142437439120778e-05
         check2 = 5.222222222222222
         assert count.mean() == check
         assert count2.mean() == check
         assert unweigthed.mean() == check2
         with pytest.raises(ValueError):
-            count = mm.blocks_count(
+            count = mm.BlocksCount(
                 self.df_tessellation, "bID", sw, "uID", weighted="yes"
             )
 
-    def test_reached(self):
-        count = mm.reached(self.df_streets, self.df_buildings, "nID", "nID")
-        area = mm.reached(
+    def test_Reached(self):
+        count = mm.Reached(self.df_streets, self.df_buildings, "nID", "nID").r
+        area = mm.Reached(
             self.df_streets,
             self.df_buildings,
             self.df_streets.nID,
             self.df_buildings.nID,
             mode="sum",
-        )
-        mean = mm.reached(self.df_streets, self.df_buildings, "nID", "nID", mode="mean")
-        std = mm.reached(self.df_streets, self.df_buildings, "nID", "nID", mode="std")
-        area_v = mm.reached(
+        ).r
+        mean = mm.Reached(
+            self.df_streets, self.df_buildings, "nID", "nID", mode="mean"
+        ).r
+        std = mm.Reached(self.df_streets, self.df_buildings, "nID", "nID", mode="std").r
+        area_v = mm.Reached(
             self.df_streets,
             self.df_buildings,
             "nID",
             "nID",
             mode="sum",
             values="fl_area",
-        )
-        mean_v = mm.reached(
+        ).r
+        mean_v = mm.Reached(
             self.df_streets,
             self.df_buildings,
             "nID",
             "nID",
             mode="mean",
             values="fl_area",
-        )
-        std_v = mm.reached(
+        ).r
+        std_v = mm.Reached(
             self.df_streets,
             self.df_buildings,
             "nID",
             "nID",
             mode="std",
             values="fl_area",
-        )
+        ).r
         sw = mm.sw_high(k=2, gdf=self.df_streets)
-        count_sw = mm.reached(self.df_streets, self.df_buildings, "nID", "nID", sw)
+        count_sw = mm.Reached(self.df_streets, self.df_buildings, "nID", "nID", sw).r
         assert max(count) == 18
         assert max(area) == 18085.45897711331
         assert max(count_sw) == 138
@@ -159,32 +157,32 @@ class TestIntensity:
         assert max(mean_v) == 7916.931385861784
         assert max(std_v) == 8995.18003493457
 
-    def test_node_density(self):
+    def test_NodeDensity(self):
         nx = mm.gdf_to_nx(self.df_streets)
         nx = mm.node_degree(nx)
         nodes, edges, W = mm.nx_to_gdf(nx, spatial_weights=True)
         sw = mm.sw_high(k=3, weights=W)
-        density = mm.node_density(nodes, edges, sw)
-        weighted = mm.node_density(
+        density = mm.NodeDensity(nodes, edges, sw).nd
+        weighted = mm.NodeDensity(
             nodes, edges, sw, weighted=True, node_degree="degree"
-        )
-        array = mm.node_density(nodes, edges, W)
+        ).nd
+        array = mm.NodeDensity(nodes, edges, W).nd
         assert density.mean() == 0.012690163074599968
         assert weighted.mean() == 0.023207675994368446
         assert array.mean() == 0.008554067995928158
 
-    def test_density(self):
+    def test_Density(self):
         sw = mm.sw_high(k=3, gdf=self.df_tessellation, ids="uID")
-        dens = mm.density(
+        dens = mm.Density(
             self.df_tessellation,
             self.df_buildings["fl_area"],
             sw,
             "uID",
             self.df_tessellation.area,
-        )
-        dens2 = mm.density(
+        ).d
+        dens2 = mm.Density(
             self.df_tessellation, self.df_buildings["fl_area"], sw, "uID"
-        )
+        ).d
         check = 1.661587
         assert dens.mean() == approx(check)
         assert dens2.mean() == approx(check)
diff --git a/tests/test_shape.py b/tests/test_shape.py
index cf9d143..0e372ea 100644
--- a/tests/test_shape.py
+++ b/tests/test_shape.py
@@ -15,76 +15,68 @@ class TestShape:
         self.df_streets = gpd.read_file(test_file_path, layer="streets")
         self.df_tessellation = gpd.read_file(test_file_path, layer="tessellation")
         self.df_buildings["height"] = np.linspace(10.0, 30.0, 144)
-        self.df_buildings["volume"] = mm.volume(self.df_buildings, "height")
+        self.df_buildings["volume"] = mm.Volume(self.df_buildings, "height").volume
 
-    def test_form_factor(self):
-        self.df_buildings["ff"] = mm.form_factor(self.df_buildings, "volume")
+    def test_FormFactor(self):
+        self.df_buildings["ff"] = mm.FormFactor(self.df_buildings, "volume").ff
         check = (self.df_buildings.geometry[0].area) / (
             self.df_buildings.volume[0] ** (2 / 3)
         )
         assert self.df_buildings["ff"][0] == check
 
-    def test_form_factor_array(self):
-        self.df_buildings["ff"] = mm.form_factor(
+        self.df_buildings["ff"] = mm.FormFactor(
             self.df_buildings,
-            mm.volume(self.df_buildings, "height"),
+            mm.Volume(self.df_buildings, "height").volume,
             areas=self.df_buildings.geometry.area,
-        )
+        ).ff
         check = (self.df_buildings.geometry[0].area) / (
             self.df_buildings.volume[0] ** (2 / 3)
         )
         assert self.df_buildings["ff"][0] == check
 
-    def test_fractal_dimension(self):
-        self.df_buildings["fd"] = mm.fractal_dimension(self.df_buildings)
+    def test_FractalDimension(self):
+        self.df_buildings["fd"] = mm.FractalDimension(self.df_buildings).fd
         check = math.log(self.df_buildings.geometry[0].length / 4) / math.log(
             self.df_buildings.geometry[0].area
         )
         assert self.df_buildings["fd"][0] == check
 
-    def test_fractal_dimension_array(self):
-        self.df_buildings["fd"] = mm.fractal_dimension(
+        self.df_buildings["fd2"] = mm.FractalDimension(
             self.df_buildings,
             areas=self.df_buildings.geometry.area,
             perimeters=self.df_buildings.geometry.length,
-        )
-        check = math.log(self.df_buildings.geometry[0].length / 4) / math.log(
-            self.df_buildings.geometry[0].area
-        )
-        assert self.df_buildings["fd"][0] == check
+        ).fd
+        assert self.df_buildings["fd2"][0] == check
 
-    def test_volume_facade_ratio(self):
+    def test_VolumeFacadeRatio(self):
         self.df_buildings["peri"] = self.df_buildings.geometry.length
-        self.df_buildings["vfr"] = mm.volume_facade_ratio(
+        self.df_buildings["vfr"] = mm.VolumeFacadeRatio(
             self.df_buildings, "height", "volume", "peri"
-        )
+        ).vfr
         check = self.df_buildings.volume[0] / (
             self.df_buildings.peri[0] * self.df_buildings.height[0]
         )
         assert self.df_buildings["vfr"][0] == check
 
-    def test_volume_facade_ratio_array(self):
         peri = self.df_buildings.geometry.length
-        volume = mm.volume(self.df_buildings, "height")
-        self.df_buildings["vfr"] = mm.volume_facade_ratio(
+        volume = mm.Volume(self.df_buildings, "height").volume
+        self.df_buildings["vfr2"] = mm.VolumeFacadeRatio(
             self.df_buildings, "height", volume, peri
-        )
-        check = self.df_buildings.volume[0] / (
-            self.df_buildings.geometry[0].length * self.df_buildings.height[0]
-        )
-        assert self.df_buildings["vfr"][0] == check
+        ).vfr
+        assert self.df_buildings["vfr2"][0] == check
 
-    def test_volume_facade_ratio_nones(self):
         self.df_buildings["peri"] = self.df_buildings.geometry.length
-        self.df_buildings["vfr"] = mm.volume_facade_ratio(self.df_buildings, "height")
-        check = self.df_buildings.volume[0] / (
-            self.df_buildings.peri[0] * self.df_buildings.height[0]
-        )
-        assert self.df_buildings["vfr"][0] == check
+        self.df_buildings["vfr3"] = mm.VolumeFacadeRatio(
+            self.df_buildings, "height"
+        ).vfr
 
-    def test_circular_compactness(self):
+        assert self.df_buildings["vfr3"][0] == check
+
+    def test_CircularCompactness(self):
         self.df_buildings["area"] = self.df_buildings.geometry.area
-        self.df_buildings["circom"] = mm.circular_compactness(self.df_buildings, "area")
+        self.df_buildings["circom"] = mm.CircularCompactness(
+            self.df_buildings, "area"
+        ).cc
         check = self.df_buildings.area[0] / (
             _circle_area(
                 list(self.df_buildings.geometry[0].convex_hull.exterior.coords)
@@ -92,78 +84,60 @@ class TestShape:
         )
         assert self.df_buildings["circom"][0] == check
 
-    def test_circular_compactness_array(self):
         area = self.df_buildings.geometry.area
-        self.df_buildings["circom"] = mm.circular_compactness(self.df_buildings, area)
-        check = self.df_buildings.area[0] / (
-            _circle_area(
-                list(self.df_buildings.geometry[0].convex_hull.exterior.coords)
-            )
-        )
-        assert self.df_buildings["circom"][0] == check
+        self.df_buildings["circom2"] = mm.CircularCompactness(
+            self.df_buildings, area
+        ).cc
+        assert self.df_buildings["circom2"][0] == check
 
-    def test_circular_compactness_nones(self):
-        self.df_buildings["circom"] = mm.circular_compactness(self.df_buildings)
-        check = self.df_buildings.area[0] / (
-            _circle_area(
-                list(self.df_buildings.geometry[0].convex_hull.exterior.coords)
-            )
-        )
-        assert self.df_buildings["circom"][0] == check
+        self.df_buildings["circom3"] = mm.CircularCompactness(self.df_buildings).cc
+        assert self.df_buildings["circom3"][0] == check
 
-    def test_square_compactness(self):
-        self.df_buildings["sqcom"] = mm.square_compactness(self.df_buildings)
+    def test_SquareCompactness(self):
+        self.df_buildings["sqcom"] = mm.SquareCompactness(self.df_buildings).sc
         check = (
             (4 * math.sqrt(self.df_buildings.geometry.area[0]))
             / (self.df_buildings.geometry.length[0])
         ) ** 2
         assert self.df_buildings["sqcom"][0] == check
 
-    def test_square_compactness_array(self):
-        self.df_buildings["sqcom"] = mm.square_compactness(
+        self.df_buildings["sqcom2"] = mm.SquareCompactness(
             self.df_buildings,
             areas=self.df_buildings.geometry.area,
             perimeters=self.df_buildings.geometry.length,
-        )
-        check = (
-            (4 * math.sqrt(self.df_buildings.geometry.area[0]))
-            / (self.df_buildings.geometry.length[0])
-        ) ** 2
-        assert self.df_buildings["sqcom"][0] == check
+        ).sc
+
+        assert self.df_buildings["sqcom2"][0] == check
 
-    def test_convexeity(self):
-        self.df_buildings["conv"] = mm.convexeity(self.df_buildings)
+    def test_Convexeity(self):
+        self.df_buildings["conv"] = mm.Convexeity(self.df_buildings).c
         check = (
             self.df_buildings.geometry.area[0]
             / self.df_buildings.geometry.convex_hull.area[0]
         )
         assert self.df_buildings["conv"][0] == check
 
-    def test_convexeity_array(self):
-        self.df_buildings["conv"] = mm.convexeity(
+        self.df_buildings["conv2"] = mm.Convexeity(
             self.df_buildings, areas=self.df_buildings.geometry.area
-        )
-        check = (
-            self.df_buildings.geometry.area[0]
-            / self.df_buildings.geometry.convex_hull.area[0]
-        )
-        assert self.df_buildings["conv"][0] == check
+        ).c
+
+        assert self.df_buildings["conv2"][0] == check
 
-    def test_courtyard_index(self):
-        cas = self.df_buildings["cas"] = mm.courtyard_area(self.df_buildings)
-        self.df_buildings["cix"] = mm.courtyard_index(self.df_buildings, "cas")
-        self.df_buildings["cix_array"] = mm.courtyard_index(
+    def test_CourtyardIndex(self):
+        cas = self.df_buildings["cas"] = mm.CourtyardArea(self.df_buildings).ca
+        self.df_buildings["cix"] = mm.CourtyardIndex(self.df_buildings, "cas").ci
+        self.df_buildings["cix_array"] = mm.CourtyardIndex(
             self.df_buildings, cas, self.df_buildings.geometry.area
-        )
+        ).ci
         check = self.df_buildings.cas[80] / self.df_buildings.geometry.area[80]
         assert self.df_buildings["cix"][80] == check
         assert self.df_buildings["cix_array"][80] == check
 
-    def test_rectangularity(self):
-        self.df_buildings["rect"] = mm.rectangularity(self.df_buildings)
-        self.df_buildings["rect_array"] = mm.rectangularity(
+    def test_Rectangularity(self):
+        self.df_buildings["rect"] = mm.Rectangularity(self.df_buildings).r
+        self.df_buildings["rect_array"] = mm.Rectangularity(
             self.df_buildings, self.df_buildings.geometry.area
-        )
+        ).r
         check = (
             self.df_buildings.geometry[0].area
             / self.df_buildings.geometry[0].minimum_rotated_rectangle.area
@@ -171,69 +145,69 @@ class TestShape:
         assert self.df_buildings["rect"][0] == check
         assert self.df_buildings["rect_array"][0] == check
 
-    def test_shape_index(self):
-        la = self.df_buildings["la"] = mm.longest_axis_length(self.df_buildings)
-        self.df_buildings["shape_index"] = mm.shape_index(self.df_buildings, "la")
-        self.df_buildings["shape_index_array"] = mm.shape_index(
+    def test_ShapeIndex(self):
+        la = self.df_buildings["la"] = mm.LongestAxisLength(self.df_buildings).lal
+        self.df_buildings["shape_index"] = mm.ShapeIndex(self.df_buildings, "la").si
+        self.df_buildings["shape_index_array"] = mm.ShapeIndex(
             self.df_buildings, la, self.df_buildings.geometry.area
-        )
+        ).si
         check = math.sqrt(self.df_buildings.area[0] / math.pi) / (
             0.5 * self.df_buildings.la[0]
         )
         assert self.df_buildings["shape_index"][0] == check
         assert self.df_buildings["shape_index_array"][0] == check
 
-    def test_corners(self):
-        self.df_buildings["corners"] = mm.corners(self.df_buildings)
+    def test_Corners(self):
+        self.df_buildings["corners"] = mm.Corners(self.df_buildings).c
         check = 24
         assert self.df_buildings["corners"][0] == check
 
-    def test_squareness(self):
-        self.df_buildings["squ"] = mm.squareness(self.df_buildings)
+    def test_Squareness(self):
+        self.df_buildings["squ"] = mm.Squareness(self.df_buildings).s
         check = round(3.7075816043359855, 4)
         assert round(self.df_buildings["squ"][0], 4) == check
 
-    def test_equivalent_rectangular_index(self):
-        self.df_buildings["eri"] = mm.equivalent_rectangular_index(self.df_buildings)
-        self.df_buildings["eri_array"] = mm.equivalent_rectangular_index(
+    def test_EquivalentRectangularIndex(self):
+        self.df_buildings["eri"] = mm.EquivalentRectangularIndex(self.df_buildings).eri
+        self.df_buildings["eri_array"] = mm.EquivalentRectangularIndex(
             self.df_buildings,
             areas=self.df_buildings.geometry.area,
             perimeters=self.df_buildings.geometry.length,
-        )
+        ).eri
         check = 0.7879229963118455
         assert self.df_buildings["eri"][0] == check
         assert self.df_buildings["eri_array"][0] == check
 
-    def test_elongation(self):
-        self.df_buildings["elo"] = mm.elongation(self.df_buildings)
+    def test_Elongation(self):
+        self.df_buildings["elo"] = mm.Elongation(self.df_buildings).e
         check = 0.9082437463675544
         assert self.df_buildings["elo"][0] == check
 
-    def test_centroid_corners(self):
-        means, devs = mm.centroid_corners(self.df_buildings)
-        self.df_buildings["ccd"] = means
-        self.df_buildings["ccddev"] = devs
+    def test_CentroidCorners(self):
+        cc = mm.CentroidCorners(self.df_buildings)
+        self.df_buildings["ccd"] = cc.mean
+        self.df_buildings["ccddev"] = cc.std
         check = 15.961531913184833
         check_devs = 3.0810634305400177
         assert self.df_buildings["ccd"][0] == check
         assert self.df_buildings["ccddev"][0] == check_devs
 
-    def test_linearity(self):
-        self.df_streets["lin"] = mm.linearity(self.df_streets)
+    def test_Linearity(self):
+        self.df_streets["lin"] = mm.Linearity(self.df_streets).linearity
         euclidean = Point(self.df_streets.geometry[0].coords[0]).distance(
             Point(self.df_streets.geometry[0].coords[-1])
         )
         check = euclidean / self.df_streets.geometry[0].length
         assert self.df_streets["lin"][0] == check
 
-    def test_compactness_weighted_axis(self):
-        self.df_buildings["cwa"] = mm.compactness_weighted_axis(self.df_buildings)
-        self.df_buildings["cwa_array"] = mm.compactness_weighted_axis(
+    def test_CompactnessWeightedAxis(self):
+        self.df_buildings["cwa"] = mm.CompactnessWeightedAxis(self.df_buildings).cwa
+        self.df_buildings["cwa_array"] = mm.CompactnessWeightedAxis(
             self.df_buildings,
             areas=self.df_buildings.geometry.area,
             perimeters=self.df_buildings.geometry.length,
-            longest_axis=mm.longest_axis_length(self.df_buildings),
-        )
+            longest_axis=mm.LongestAxisLength(self.df_buildings).lal,
+        ).cwa
         check = 26.32772969906327
         assert self.df_buildings["cwa"][0] == check
         assert self.df_buildings["cwa_array"][0] == check
diff --git a/tests/test_utils.py b/tests/test_utils.py
index eeb699e..cc0e1d8 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -73,9 +73,9 @@ class TestUtils:
             mm.nx_to_gdf(G)
 
     def test_limit_range(self):
-        assert mm.limit_range(range(10), rng=(25, 75)) == [2, 3, 4, 5, 6, 7]
-        assert mm.limit_range(range(10), rng=(10, 90)) == [1, 2, 3, 4, 5, 6, 7, 8]
-        assert mm.limit_range([0, 1], rng=(25, 75)) == [0, 1]
+        assert list(mm.limit_range(range(10), rng=(25, 75))) == [2, 3, 4, 5, 6, 7]
+        assert list(mm.limit_range(range(10), rng=(10, 90))) == [1, 2, 3, 4, 5, 6, 7, 8]
+        assert list(mm.limit_range([0, 1], rng=(25, 75))) == [0, 1]
 
     def test_preprocess(self):
         test_file_path2 = mm.datasets.get_path("tests")