Repository: martinfleis/momepy
Pull request #12 (2019-06-03 14:43:31)

ENH: Added diversity characters, wall, covered area, wall (plus minor refactor)
Added couple of functions. At the same time minor terminology change and refactor of limit_range(). Plus changes on travis. (Divide into more PR next time.)

codecov[bot]: # [Codecov](https://codecov.io/gh/martinfleis/momepy/pull/12?src=pr&el=h1) Report
> Merging [#12](https://codecov.io/gh/martinfleis/momepy/pull/12?src=pr&el=desc) into [master](https://codecov.io/gh/martinfleis/momepy/commit/84a05b764f840317a8740a9de406f6333c62f27e?src=pr&el=desc) will **increase** coverage by `2.38%`.
> The diff coverage is `90.76%`.

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

```diff
@@            Coverage Diff             @@
##           master      #12      +/-   ##
==========================================
+ Coverage   86.03%   88.41%   +2.38%     
==========================================
  Files           9        9              
  Lines        2026     2107      +81     
==========================================
+ Hits         1743     1863     +120     
+ Misses        283      244      -39
```


| [Impacted Files](https://codecov.io/gh/martinfleis/momepy/pull/12?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [momepy/elements.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L2VsZW1lbnRzLnB5) | `82.35% <100%> (ø)` | :arrow_up: |
| [momepy/distribution.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L2Rpc3RyaWJ1dGlvbi5weQ==) | `96.48% <100%> (ø)` | :arrow_up: |
| [momepy/utils.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L3V0aWxzLnB5) | `87.89% <100%> (-0.38%)` | :arrow_down: |
| [momepy/intensity.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L2ludGVuc2l0eS5weQ==) | `93.1% <100%> (+13.96%)` | :arrow_up: |
| [momepy/diversity.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L2RpdmVyc2l0eS5weQ==) | `84.31% <83.5%> (+66.13%)` | :arrow_up: |
| [momepy/dimension.py](https://codecov.io/gh/martinfleis/momepy/pull/12/diff?src=pr&el=tree#diff-bW9tZXB5L2RpbWVuc2lvbi5weQ==) | `97.65% <96.29%> (-0.33%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/martinfleis/momepy/pull/12?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/12?src=pr&el=footer). Last update [84a05b7...0b12eee](https://codecov.io/gh/martinfleis/momepy/pull/12?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/.travis.yml b/.travis.yml
index 98da078..eb3f862 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,20 +2,17 @@ language: python
 
 sudo: false
 
+dist: xenial   # required for Python >= 3.7
+
 python:
-  - "2.7"
-  - "3.5"
   - "3.6"
+  - "3.7"
 
 install:
   - sudo apt-get update
   # We do this conditionally because it saves us some downloading if the
   # version is the same.
-  - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then
-      wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh;
-    else
-      wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
-    fi
+  - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
   - bash miniconda.sh -b -p $HOME/miniconda
   - export PATH="$HOME/miniconda/bin:$PATH"
   - hash -r
@@ -28,7 +25,7 @@ install:
 
   - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION
   - source activate test-environment
-  - conda install geopandas libpysal osmnx tqdm pytest pytest pytest-cov codecov
+  - conda install geopandas libpysal osmnx tqdm pysal mapclassify pytest pytest-cov codecov
   - python setup.py install
   - conda list
 
diff --git a/docs/api.rst b/docs/api.rst
index 10023c1..0a4ffd7 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -26,12 +26,14 @@ dimension
 
    area
    courtyard_area
+   covered_area
    floor_area
    longest_axis_length
    mean_character
    perimeter
    street_profile
    volume
+   wall
    weighted_character
 
 shape
@@ -86,13 +88,24 @@ intensity
    frequency
    gross_density
 
+diversity
+---------
+.. autosummary::
+  :toctree: generated/
+
+  rng
+  simpson
+  theil
+
 utilities
 ---------
 .. autosummary::
    :toctree: generated/
 
-   Queen_higher
    gdf_to_nx
+   limit_range
+   multi2single
    nx_to_gdf
+   preprocess
    unique_id
-   from_os
+   Queen_higher
diff --git a/momepy/dimension.py b/momepy/dimension.py
index 21cd38e..2d3ef41 100644
--- a/momepy/dimension.py
+++ b/momepy/dimension.py
@@ -299,7 +299,7 @@ def longest_axis_length(objects):
     return series
 
 
-def mean_character(objects, spatial_weights=None, values=None, order=3, mode=None):
+def mean_character(objects, spatial_weights=None, values=None, order=3, rng=None):
     """
     Calculates the mean of a character within k steps of morphological tessellation
 
@@ -319,8 +319,9 @@ def mean_character(objects, spatial_weights=None, values=None, order=3, mode=Non
         the name of the dataframe column, np.array, or pd.Series where is stored character value. If set to None, function will use area.
     order : int
         order of Queen contiguity
-    mode : str
-        mode of calculation. None will use all values, `iq` will use values within interquartile range, `id` will use values within interdecile range
+    rng : Two-element sequence containing floats in range of [0,100], optional
+        Percentiles over which to compute the range. Each must be
+        between 0 and 100, inclusive. The order of the elements is not important.
 
     Returns
     -------
@@ -384,9 +385,9 @@ def mean_character(objects, spatial_weights=None, values=None, order=3, mode=Non
                 values_list.append(row.geometry.area)
             else:
                 values_list = [row.geometry.area]
-        if mode:
+        if rng:
             from momepy import limit_range
-            values_list = limit_range(values_list, mode=mode)
+            values_list = limit_range(values_list, rng=rng)
         results_list.append(sum(values_list) / len(values_list))
 
     series = pd.Series(results_list)
@@ -693,12 +694,127 @@ def weighted_character(objects, tessellation, characters, unique_id, spatial_wei
 
     print('Weighted {} calculated.'.format(characters))
     return series
-# to be deleted, keep at the end
 
-# path = "/Users/martin/Dropbox/StrathUni/PhD/Papers/Voronoi tesselation/Data/Zurich/Final data/Voronoi/test/voronoi_10.shp"
-# objects = gpd.read_file(path)
 
-# longest_axis_length2(objects='longest_axis')
+def covered_area(objects, spatial_weights=None):
+    """
+    Calculates the area covered by k steps of morphological tessellation
+
+    Total area covered within k topological steps defined in spatial_weights.
+
+    .. math::
+
+
+    Parameters
+    ----------
+    objects : GeoDataFrame
+        GeoDataFrame containing morphological tessellation
+    spatial_weights : libpysal.weights, optional
+        spatial weights matrix - If None, Queen contiguity matrix will be calculated
+        based on objects.
+
+    Returns
+    -------
+    Series
+        Series containing resulting values.
+
+    References
+    ----------
+
+    Examples
+    --------
+
+    """
+    # define empty list for results
+    results_list = []
+
+    print('Calculating covered area...')
+
+    if spatial_weights is None:
+        print('Generating weights matrix (Queen)...')
+        from libpysal.weights import Queen
+        # matrix to define area of analysis (more steps)
+        spatial_weights = Queen.from_dataframe(objects, silence_warnings=True)
+    else:
+        if not all(objects.index == range(len(objects))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
+
+    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
+        neighbours = spatial_weights.neighbors[index]
+        neighbours.append(index)
+        areas = objects.iloc[neighbours].geometry.area
+        results_list.append(sum(areas))
+
+    series = pd.Series(results_list)
+
+    print('Covered area calculated.')
+    return series
+
+
+def wall(objects, spatial_weights=None):
+    """
+    Calculate the perimeter wall length the joined structure.
+
+    Parameters
+    ----------
+    objects : GeoDataFrame
+        GeoDataFrame containing objects to analyse
+    spatial_weights : libpysal.weights, optional
+        spatial weights matrix - If None, Queen contiguity matrix will be calculated
+        based on objects. It is to denote adjacent buildings.
+
+    Returns
+    -------
+    Series
+        Series containing resulting values.
+
+    Notes
+    -----
+    Script is not optimised at all, so it is currently extremely slow.
+    """
+    # define empty list for results
+    results_list = []
+
+    print('Calculating perimeter wall length...')
+
+    if not all(objects.index == range(len(objects))):
+        raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
-# objects.head
-# objects.to_file(path)
+    # 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(objects, silence_warnings=True)
+        print('Spatial weights ready...')
+
+    # dict to store walls for each uID
+    walls = {}
+
+    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
+        # if the id is already present in walls, continue (avoid repetition)
+        if index in walls:
+            continue
+        else:
+            to_join = [index]  # list of indices which should be joined together
+            neighbours = []  # list of neighbours
+            weights = spatial_weights.neighbors[index]  # neighbours from spatial weights
+            for w in weights:
+                neighbours.append(w)  # make a list from weigths
+
+            for n in neighbours:
+                while n not in to_join:  # until there is some neighbour which is not in to_join
+                    to_join.append(n)
+                    weights = spatial_weights.neighbors[n]
+                    for w in weights:
+                        neighbours.append(w)  # extend neighbours by neighbours of neighbours :)
+            joined = objects.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(objects.iterrows(), total=objects.shape[0]):
+        results_list.append(walls[index])
+
+    series = pd.Series(results_list)
+    print('Perimeter wall length calculated.')
+    return series
diff --git a/momepy/distribution.py b/momepy/distribution.py
index c981917..e2e284f 100644
--- a/momepy/distribution.py
+++ b/momepy/distribution.py
@@ -323,7 +323,7 @@ def cell_alignment(objects, tessellation, orientation_column, cell_orientation_c
     return series
 
 
-def alignment(objects, orientation_column, tessellation, unique_id, weights_matrix=None):
+def alignment(objects, orientation_column, tessellation, unique_id, spatial_weights=None):
     """
     Calculate the mean deviation of solar orientation of objects on adjacent cells from object
 
@@ -337,12 +337,12 @@ def alignment(objects, orientation_column, tessellation, unique_id, weights_matr
     orientation_column : str
         name of the column where is stored object orientation value
     tessellation : GeoDataFrame
-        GeoDataFrame containing morphological tessellation - source of weights_matrix.
+        GeoDataFrame containing morphological tessellation - source of spatial_weights.
         It is crucial to use exactly same input as was used durign the calculation of weights matrix.
-        If weights_matrix is None, tessellation is used to calulate it.
+        If spatial_weights is None, tessellation is used to calulate it.
     unique_id : str
         name of the column with unique id
-    weights_matrix : libpysal.weights, optional
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on tessellation
 
@@ -359,16 +359,16 @@ def alignment(objects, orientation_column, tessellation, unique_id, weights_matr
     if not all(tessellation.index == range(len(tessellation))):
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Calculating spatial weights...')
         from libpysal.weights import Queen
-        weights_matrix = Queen.from_dataframe(tessellation)
+        spatial_weights = Queen.from_dataframe(tessellation)
         print('Spatial weights ready...')
 
     # iterating over rows one by one
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
         uid = tessellation.loc[tessellation[unique_id] == row[unique_id]].index[0]
-        neighbours = weights_matrix.neighbors[uid]
+        neighbours = spatial_weights.neighbors[uid]
         neighbours_ids = []
 
         for n in neighbours:
@@ -396,7 +396,7 @@ def alignment(objects, orientation_column, tessellation, unique_id, weights_matr
     return series
 
 
-def neighbour_distance(objects, tessellation, unique_id, weights_matrix=None):
+def neighbour_distance(objects, tessellation, unique_id, spatial_weights=None):
     """
     Calculate the mean distance to buildings on adjacent cells
 
@@ -408,12 +408,12 @@ def neighbour_distance(objects, tessellation, unique_id, weights_matrix=None):
     objects : GeoDataFrame
         GeoDataFrame containing objects to analyse
     tessellation : GeoDataFrame
-        GeoDataFrame containing morphological tessellation - source of weights_matrix.
+        GeoDataFrame containing morphological tessellation - source of spatial_weights.
         It is crucial to use exactly same input as was used durign the calculation of weights matrix.
-        If weights_matrix is None, tessellation is used to calulate it.
+        If spatial_weights is None, tessellation is used to calulate it.
     unique_id : str
         name of the column with unique id
-    weights_matrix : libpysal.weights, optional
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on tessellation
 
@@ -435,16 +435,16 @@ def neighbour_distance(objects, tessellation, unique_id, weights_matrix=None):
     if not all(tessellation.index == range(len(tessellation))):
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Calculating spatial weights...')
         from libpysal.weights import Queen
-        weights_matrix = Queen.from_dataframe(tessellation)
+        spatial_weights = Queen.from_dataframe(tessellation)
         print('Spatial weights ready...')
 
     # iterating over rows one by one
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
         uid = tessellation.loc[tessellation[unique_id] == row[unique_id]].index[0]
-        neighbours = weights_matrix.neighbors[uid]
+        neighbours = spatial_weights.neighbors[uid]
 
         neighbours_ids = tessellation.iloc[neighbours][unique_id]
         building_neighbours = objects.loc[objects[unique_id].isin(neighbours_ids)]
@@ -459,11 +459,11 @@ def neighbour_distance(objects, tessellation, unique_id, weights_matrix=None):
     return series
 
 
-def mean_interbuilding_distance(objects, tessellation, unique_id, weights_matrix=None, weights_matrix_higher=None, order=3):
+def mean_interbuilding_distance(objects, tessellation, unique_id, spatial_weights=None, spatial_weights_higher=None, order=3):
     """
     Calculate the mean interbuilding distance within x topological steps
 
-    Interbuilding distances are calculated between buildings on adjacent cells based on `weights_matrix`.
+    Interbuilding distances are calculated between buildings on adjacent cells based on `spatial_weights`.
 
     .. math::
 
@@ -473,15 +473,15 @@ def mean_interbuilding_distance(objects, tessellation, unique_id, weights_matrix
     objects : GeoDataFrame
         GeoDataFrame containing objects to analyse
     tessellation : GeoDataFrame
-        GeoDataFrame containing morphological tessellation - source of weights_matrix and weights_matrix_higher.
-        It is crucial to use exactly same input as was used durign the calculation of weights matrix and weights_matrix_higher.
-        If weights_matrix or weights_matrix_higher is None, tessellation is used to calulate it.
+        GeoDataFrame containing morphological tessellation - source of spatial_weights and spatial_weights_higher.
+        It is crucial to use exactly same input as was used durign the calculation of weights matrix and spatial_weights_higher.
+        If spatial_weights or spatial_weights_higher is None, tessellation is used to calulate it.
     unique_id : str
         name of the column with unique id
-    weights_matrix : libpysal.weights, optional
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on tessellation
-    weights_matrix_higher : libpysal.weights, optional
+    spatial_weights_higher : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity of higher order will be calculated
         based on tessellation
     order : int
@@ -500,24 +500,24 @@ def mean_interbuilding_distance(objects, tessellation, unique_id, weights_matrix
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
     print('Calculating mean interbuilding distances...')
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Generating weights matrix (Queen)...')
         from libpysal.weights import Queen
         # matrix to capture interbuilding relationship
-        weights_matrix = Queen.from_dataframe(tessellation)
+        spatial_weights = Queen.from_dataframe(tessellation)
 
-    if weights_matrix_higher is None:
+    if spatial_weights_higher is None:
         print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
         from momepy import Queen_higher
         # matrix to define area of analysis (more steps)
-        weights_matrix_higher = Queen_higher(k=order, geodataframe=tessellation)
+        spatial_weights_higher = Queen_higher(k=order, geodataframe=tessellation)
 
     # define empty list for results
     results_list = []
 
     print('Generating adjacency matrix based on weights matrix...')
     # define adjacency list from lipysal
-    adj_list = weights_matrix.to_adjlist()
+    adj_list = spatial_weights.to_adjlist()
     adj_list['distance'] = -1
 
     print('Computing interbuilding distances...')
@@ -540,7 +540,7 @@ def mean_interbuilding_distance(objects, tessellation, unique_id, weights_matrix
         # id to match spatial weights
         uid = tessellation.loc[tessellation[unique_id] == row[unique_id]].index[0]
         # define neighbours based on weights matrix defining analysis area
-        neighbours = weights_matrix_higher.neighbors[uid]
+        neighbours = spatial_weights_higher.neighbors[uid]
         neighbours.append(uid)
         if neighbours:
             selection = adj_list[adj_list.focal.isin(neighbours)][adj_list.neighbor.isin(neighbours)]
@@ -638,7 +638,7 @@ def neighbouring_street_orientation_deviation(objects):
     return series
 
 
-def building_adjacency(objects, tessellation, weights_matrix=None, weights_matrix_higher=None, order=3, unique_id='uID'):
+def building_adjacency(objects, tessellation, spatial_weights=None, spatial_weights_higher=None, order=3, unique_id='uID'):
     """
     Calculate the level of building adjacency
 
@@ -653,13 +653,13 @@ def building_adjacency(objects, tessellation, weights_matrix=None, weights_matri
     objects : GeoDataFrame
         GeoDataFrame containing objects to analyse
     tessellation : GeoDataFrame
-        GeoDataFrame containing morphological tessellation - source of weights_matrix and weights_matrix_higher.
-        It is crucial to use exactly same input as was used durign the calculation of weights matrix and weights_matrix_higher.
-        If weights_matrix or weights_matrix_higher is None, tessellation is used to calulate it.
-    weights_matrix : libpysal.weights, optional
+        GeoDataFrame containing morphological tessellation - source of spatial_weights and spatial_weights_higher.
+        It is crucial to use exactly same input as was used durign the calculation of weights matrix and spatial_weights_higher.
+        If spatial_weights or spatial_weights_higher is None, tessellation is used to calulate it.
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on tessellation
-    weights_matrix_higher : libpysal.weights, optional
+    spatial_weights_higher : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity of higher order will be calculated
         based on tessellation
     order : int
@@ -684,17 +684,17 @@ def building_adjacency(objects, tessellation, weights_matrix=None, weights_matri
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
     # if weights matrix is not passed, generate it from objects
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Calculating spatial weights...')
         from libpysal.weights import Queen
-        weights_matrix = Queen.from_dataframe(objects, silence_warnings=True)
+        spatial_weights = Queen.from_dataframe(objects, silence_warnings=True)
         print('Spatial weights ready...')
 
-    if weights_matrix_higher is None:
+    if spatial_weights_higher is None:
         print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
         from momepy import Queen_higher
         # matrix to define area of analysis (more steps)
-        weights_matrix_higher = Queen_higher(k=order, geodataframe=tessellation)
+        spatial_weights_higher = Queen_higher(k=order, geodataframe=tessellation)
 
     print('Generating dictionary of built-up patches...')
     # dict to store nr of courtyards for each uID
@@ -708,14 +708,14 @@ def building_adjacency(objects, tessellation, weights_matrix=None, weights_matri
         else:
             to_join = [index]  # list of indices which should be joined together
             neighbours = []  # list of neighbours
-            weights = weights_matrix.neighbors[index]  # neighbours from spatial weights
+            weights = spatial_weights.neighbors[index]  # neighbours from spatial weights
             for w in weights:
                 neighbours.append(w)  # make a list from weigths
 
             for n in neighbours:
                 while n not in to_join:  # until there is some neighbour which is not in to_join
                     to_join.append(n)
-                    weights = weights_matrix.neighbors[n]
+                    weights = spatial_weights.neighbors[n]
                     for w in weights:
                         neighbours.append(w)  # extend neighbours by neighbours of neighbours :)
             for b in to_join:
@@ -725,7 +725,7 @@ def building_adjacency(objects, tessellation, weights_matrix=None, weights_matri
     print('Calculating adjacency within k steps...')
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
         uid = tessellation.loc[tessellation[unique_id] == row[unique_id]].index[0]
-        neighbours = weights_matrix_higher.neighbors[uid]
+        neighbours = spatial_weights_higher.neighbors[uid]
 
         neighbours_ids = tessellation.iloc[neighbours][unique_id]
         neighbours_ids = neighbours_ids.append(pd.Series(row[unique_id], index=[index]))
@@ -742,7 +742,7 @@ def building_adjacency(objects, tessellation, weights_matrix=None, weights_matri
     return series
 
 
-def neighbours(objects, weights_matrix=None, weighted=False):
+def neighbours(objects, spatial_weights=None, weighted=False):
     """
     Calculate the number of topological neighbours of each object.
 
@@ -756,7 +756,7 @@ def neighbours(objects, weights_matrix=None, weighted=False):
     ----------
     objects : GeoDataFrame
         GeoDataFrame containing objects to analyse
-    weights_matrix : libpysal.weights (default None)
+    spatial_weights : libpysal.weights (default None)
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on tessellation
     weighted : bool (default False)
@@ -778,19 +778,19 @@ def neighbours(objects, weights_matrix=None, weighted=False):
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
     # if weights matrix is not passed, generate it from objects
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Calculating spatial weights...')
         from libpysal.weights import Queen
-        weights_matrix = Queen.from_dataframe(objects, silence_warnings=True)
+        spatial_weights = Queen.from_dataframe(objects, silence_warnings=True)
         print('Spatial weights ready...')
 
     print('Calculating neighbours...')
     neighbours = []
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
         if weighted is True:
-            neighbours.append(len(weights_matrix.neighbors[index]) / row.geometry.length)
+            neighbours.append(spatial_weights.cardinalities[index] / row.geometry.length)
         else:
-            neighbours.append(len(weights_matrix.neighbors[index]))
+            neighbours.append(spatial_weights.cardinalities[index])
 
     series = pd.Series(neighbours)
     print('Neighbours calculated.')
diff --git a/momepy/diversity.py b/momepy/diversity.py
index 2c99916..5cc1460 100644
--- a/momepy/diversity.py
+++ b/momepy/diversity.py
@@ -5,189 +5,282 @@
 # definitons of diversity characters
 
 from tqdm import tqdm  # progress bar
-import numpy as np
-from .intensity import radius
 import pandas as pd
-# import powerlaw
+import scipy as sp
 
 
-'''
-Gini index calculation. Gini script by https://github.com/oliviaguest/gini.
+def rng(objects, values, spatial_weights=None, order=None, rng=(0, 100), **kwargs):
+    """
+    Calculates the range of values within k steps of morphological tessellation
 
-There should be docs.
-'''
+    Uses `scipy.stats.iqr` under the hood.
 
+    .. math::
 
-def gini(array):
-    """Calculate the Gini coefficient of a numpy array."""
-    # based on bottom eq:
-    # http://www.statsdirect.com/help/generatedimages/equations/equation154.svg
-    # from:
-    # http://www.statsdirect.com/help/default.htm#nonparametric_methods/gini.htm
-    # All values are treated equally, arrays must be 1d:
-    array = array.flatten()
-    if np.amin(array) < 0:
-        # Values cannot be negative:
-        array -= np.amin(array)
-    # Values cannot be 0:
-    array += 0.0000001
-    # Values must be sorted:
-    array = np.sort(array)
-    # Index per array element:
-    index = np.arange(1, array.shape[0] + 1)
-    # Number of array elements:
-    n = array.shape[0]
-    # Gini coefficient:
-    return ((np.sum((2 * index - n - 1) * array)) / (n * np.sum(array)))
 
+    Parameters
+    ----------
+    objects : GeoDataFrame
+        GeoDataFrame containing morphological tessellation
+    values : str, list, np.array, pd.Series
+        the name of the dataframe column, np.array, or pd.Series where is stored character value.
+    spatial_weights : libpysal.weights, optional
+        spatial weights matrix - If None, Queen contiguity matrix of set order will be calculated
+        based on objects.
+    order : int
+        order of Queen contiguity
+    rng : Two-element sequence containing floats in range of [0,100], optional
+        Percentiles over which to compute the range. Each must be
+        between 0 and 100, inclusive. The order of the elements is not important.
+    **kwargs : keyword arguments
+        optional arguments for `scipy.stats.iqr`
 
-def gini_index(objects, look_for, source, column_name, id_column='uID'):
+    Returns
+    -------
+    Series
+        Series containing resulting values.
 
-    print('Calculating Gini index...')
+    References
+    ----------
 
-    objects_centroids = objects.copy()
-    objects_centroids['geometry'] = objects_centroids.centroid
+    Examples
+    --------
 
-    look_for_centroids = objects.copy()
-    look_for_centroids['geometry'] = look_for_centroids.centroid
+    """
+    # define empty list for results
+    results_list = []
+
+    print('Calculating range...')
+
+    if spatial_weights is None:
+        print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
+        from momepy import Queen_higher
+        # matrix to define area of analysis (more steps)
+        spatial_weights = Queen_higher(k=order, geodataframe=objects)
+    else:
+        if not all(objects.index == range(len(objects))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
+
+    if values is not None:
+        if not isinstance(values, str):
+            objects['mm_v'] = values
+            values = 'mm_v'
+
+    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
+        neighbours = spatial_weights.neighbors[index]
+        values_list = objects.iloc[neighbours][values].tolist()
+        if values_list:
+            values_list.append(row[values])
+        else:
+            values_list = [row[values]]
+
+        results_list.append(sp.stats.iqr(values_list, rng=rng, **kwargs))
+
+    series = pd.Series(results_list)
+
+    if 'mm_v' in objects.columns:
+        objects.drop(columns=['mm_v'], inplace=True)
+
+    print('Range calculated.')
+    return series
+
+
+def theil(objects, values, spatial_weights=None, order=None, rng=(0, 100)):
+    """
+    Calculates the Theil measure of inequality of values within k steps of morphological tessellation
+
+    Uses `pysal.explore.inequality.theil.Theil` under the hood.
+
+    .. math::
+
+
+    Parameters
+    ----------
+    objects : GeoDataFrame
+        GeoDataFrame containing morphological tessellation
+    values : str, list, np.array, pd.Series
+        the name of the dataframe column, np.array, or pd.Series where is stored character value.
+    spatial_weights : libpysal.weights, optional
+        spatial weights matrix - If None, Queen contiguity matrix of set order will be calculated
+        based on objects.
+    order : int
+        order of Queen contiguity
+    rng : Two-element sequence containing floats in range of [0,100], optional
+        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.
+
+    References
+    ----------
+
+    Examples
+    --------
+
+    """
+    from pysal.explore.inequality.theil import Theil
     # define empty list for results
     results_list = []
 
-    for index, row in tqdm(objects_centroids.iterrows(), total=objects_centroids.shape[0]):
-        neighbours = radius(look_for_centroids, row['geometry'], 400)
-        rows = objects.iloc[neighbours]
-        values = rows[[source]].values
-        results_list.append(gini(values))
+    print('Calculating Theil index...')
+
+    if spatial_weights is None:
+        print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
+        from momepy import Queen_higher
+        # matrix to define area of analysis (more steps)
+        spatial_weights = Queen_higher(k=order, geodataframe=objects)
+    else:
+        if not all(objects.index == range(len(objects))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
+
+    if values is not None:
+        if not isinstance(values, str):
+            objects['mm_v'] = values
+            values = 'mm_v'
+
+    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
+        neighbours = spatial_weights.neighbors[index]
+        values_list = objects.iloc[neighbours][values].tolist()
+        if values_list:
+            values_list.append(row[values])
+        else:
+            values_list = [row[values]]
+
+        if rng:
+            from momepy import limit_range
+            values_list = limit_range(values_list, rng=rng)
+        results_list.append(Theil(values_list).T)
 
     series = pd.Series(results_list)
 
-    print('Gini index calculated.')
+    if 'mm_v' in objects.columns:
+        objects.drop(columns=['mm_v'], inplace=True)
+
+    print('Theil index calculated.')
     return series
 
 
-'''
-Power law calculation.
-
-Docs missing.
-'''
-
-
-def power_law(objects, look_for, source, id_column='uID'):
-    # define new column
-
-    print('Calculating Power law...')
-
-    objects_centroids = objects.copy()
-    objects_centroids['geometry'] = objects_centroids.centroid
-
-    look_for_centroids = objects.copy()
-    look_for_centroids['geometry'] = look_for_centroids.centroid
-
-    for index, row in tqdm(objects_centroids.iterrows(), total=objects_centroids.shape[0]):
-        neighbours = radius(look_for_centroids, row['geometry'], 400)
-        rows = objects.iloc[neighbours]
-        values = rows[[source]]
-        values_list = values[source].tolist()
-        results = powerlaw.Fit(values_list)
-        objects.loc[index, 'pw_alpha'] = results.power_law.alpha
-        objects.loc[index, 'pw_xmin'] = results.power_law.xmin
-        # R, p = results.distribution_compare('power_law', 'lognormal')
-        # objects.loc[index, 'pw_r'] = R
-        # objects.loc[index, 'pw_p'] = p
-    return objects
-
-
-'''
-Spatial autocorrelation.
-
-Docs missing.
-'''
-# from libpysal.weights import DistanceBand
-# import esda
-#
-#
-# def moran_i_local(objects, source, column_name):
-#     print("Calculating local Moran's I...")
-#     print('Calculating weight matrix... (It might take a while.)')
-#     W = DistanceBand.from_dataframe(objects, 400)  # weight matrix
-#     y = objects[[source]]
-#     print('Calculating...')
-#     lm = esda.Moran_Local(y, W)
-#     objects[column_name] = lm.Is
-#     print("Local Moran's I calculated.")
-
-
-# import geopandas as gpd
-# import esda
-# from pysal.lib import weights, examples
-#
-#
-#
-# path = "/Users/martin/Dropbox/StrathUni/PhD/Sample Data/Prague/P7/p7_voro_single2.shp"
-# test = gpd.read_file(path)
-# test.columns
-# # testdrop = test.drop(test.index[[1998, 2070, 4216]])  # ignore islands
-# W = pysal.lib.weights.Distance.DistanceBand.from_dataframe(test, 400)  # weight matrix 'queen'
-# Wqueen
-# Wqueen = Queen.from_dataframe(test)  # weight matrix 'rook'
-# y = test[['ptcAre']]  # column to measure LISA
-# # y.drop(y.index[[1998, 2070, 4216]])
-# lm = esda.Moran_Local(y, Wqueen)  # calculate LISA Local Moran's I
-#
-# # save LISA values (need to figure out which are useful)
-# test['400local'] = lm.Is
-# test['sig'] = lm.p_sim
-# test['lm_q'] = lm.q
-# test['lm_sim'] = lm.z_sim
-#
-# gg = pysal.G_Local(y, Wrook)  # calculate LISA Gettis Ord G
-#
-# test['5local']
-# # save LISA values (need to figure out which are useful)
-# testdrop['g_Gs'] = gg.Gs
-# testdrop['g_EGs'] = gg.EGs
-# testdrop['g_Zs'] = gg.Zs
-# testdrop['g_psim'] = gg.p_sim
-#
-#
-# path2 = "/Users/martin/Strathcloud/Personal Folders/Test data/Prague/p7_voro_single2.shp"
-# test.to_file(path)
-
-
-'''
-HDBSCAN clustering. Should be in different file than diversity.
-
-Remove.
-'''
-
-# import pandas as pd
-# values = pd.read_csv('/Users/martin/Strathcloud/Personal Folders/Test data/Prague/Prague/CSV/prg.csv')
-#
-# import hdbscan
-# values.columns
-#
-# # select columns for analysis
-# clustering = values[['pdbHei', 'pdbAre', 'uID',
-#                      'pdbFlA', 'pdbVol', 'pdbPer', 'pdbCoA', 'pdbLAL', 'psbFoF', 'psbFra',
-#                      'psbVFR', 'psbCom', 'psbCom2', 'psbCon', 'psbCoI', 'psbCor', 'psbShI',
-#                      'psbSqu', 'psbERI', 'psbElo', 'ptbOri', 'pdcLAL', 'ptcAre',
-#                      'pscCom', 'pscElo', 'pivCAR', 'vicFre', 'vvcGAr', 'vvvGCA', 'uvvMCA',
-#                      'pivFAR', 'bID', 'bdkAre', 'bdkPer', 'bdkSec', 'bskElo', 'bskFra',
-#                      'bskCom', 'bskCom2', 'bskCon', 'bskShI', 'bskERI', 'btkOri']]
-# clustering[['pivFAR', 'bID', 'bdkAre', 'bdkPer', 'bdkSec', 'bskElo', 'bskFra']].isnull().sum().sum()
-# clustering.isnull().sum().sum()
-# nan_rows = clustering[clustering['bskFra'].isnull()]
-# nan_rows
-# # make clusterer
-# clusterer = hdbscan.HDBSCAN(min_cluster_size=3)
-#
-# # cluster data
-# clusterer.fit(clustering)
-# # get number of clusters
-# clusterer.labels_.max()
-# # save cluster labels to geoDataFrame
-# values['cluster'] = clusterer.labels_
-# values['cluster']
-# togis = values[['uID', 'cluster']]
-# togis.to_csv('/Users/martin/Strathcloud/Personal Folders/Test data/Prague/Prague/CSV/prg_clustered.csv')
+def simpson(objects, values, spatial_weights=None, order=None, binning='HeadTail_Breaks', **classification_kwds):
+    """
+    Calculates the Simpson\'s diversity index of values within k steps of morphological tessellation
+
+    Uses `mapclassify.classifiers` under the hood for binning.
+
+    .. math::
+
+
+    Parameters
+    ----------
+    objects : GeoDataFrame
+        GeoDataFrame containing morphological tessellation
+    values : str, list, np.array, pd.Series
+        the name of the dataframe column, np.array, or pd.Series where is stored character value.
+    spatial_weights : libpysal.weights, optional
+        spatial weights matrix - If None, Queen contiguity matrix of set order will be calculated
+        based on objects.
+    order : int
+        order of Queen contiguity
+    binning : str
+        One of mapclassify classification schemes
+        Options are Box_Plot, Equal_Interval, Fisher_Jenks,
+        Fisher_Jenks_Sampled, HeadTail_Breaks, Jenks_Caspall,
+        Jenks_Caspall_Forced, Jenks_Caspall_Sampled, Max_P_Classifier,
+        Maximum_Breaks, Natural_Breaks, Quantiles, Percentiles, Std_Mean,
+        User_Defined
+    **classification_kwds : dict
+        Keyword arguments for classification scheme
+        For details see mapclassify documentation:
+        https://mapclassify.readthedocs.io/en/latest/api.html
+
+    Returns
+    -------
+    Series
+        Series containing resulting values.
+
+    References
+    ----------
+
+    Examples
+    --------
+
+    """
+    def simpson_di(data):
+
+        """ Given a hash { 'species': count } , returns the Simpson Diversity Index
+
+        >>> simpson_di({'a': 10, 'b': 20, 'c': 30,})
+        0.3888888888888889
+
+        https://gist.github.com/martinjc/f227b447791df8c90568
+        """
+
+        def p(n, N):
+            """ Relative abundance """
+            if n == 0:
+                return 0
+            else:
+                return float(n) / N
+
+        N = sum(data.values())
+
+        return sum(p(n, N)**2 for n in data.values() if n != 0)
+
+    try:
+        import mapclassify.classifiers
+    except ImportError:
+        raise ImportError(
+            "The 'mapclassify' package is required to use the 'scheme' "
+            "keyword")
+
+    schemes = {}
+    for classifier in mapclassify.classifiers.CLASSIFIERS:
+        schemes[classifier.lower()] = getattr(mapclassify.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 = []
+
+    print('Calculating Simpson\'s diversity index...')
+
+    if spatial_weights is None:
+        print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
+        from momepy import Queen_higher
+        # matrix to define area of analysis (more steps)
+        spatial_weights = Queen_higher(k=order, geodataframe=objects)
+    else:
+        if not all(objects.index == range(len(objects))):
+            raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
+
+    if values is not None:
+        if not isinstance(values, str):
+            objects['mm_v'] = values
+            values = 'mm_v'
+
+    bins = schemes[binning](objects[values], **classification_kwds).bins
+
+    for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
+        neighbours = spatial_weights.neighbors[index]
+        values_list = objects.iloc[neighbours][values].tolist()
+        if values_list:
+            values_list.append(row[values])
+        else:
+            values_list = [row[values]]
+
+        sample_bins = mapclassify.classifiers.User_Defined(values_list, bins)
+        counts = dict(zip(bins, sample_bins.counts))
+        results_list.append(simpson_di(counts))
+
+    series = pd.Series(results_list)
+
+    if 'mm_v' in objects.columns:
+        objects.drop(columns=['mm_v'], inplace=True)
+
+    print('Simpson\'s diversity index calculated.')
+    return series
diff --git a/momepy/elements.py b/momepy/elements.py
index 6eda887..9d023b3 100644
--- a/momepy/elements.py
+++ b/momepy/elements.py
@@ -582,7 +582,7 @@ def blocks(cells, streets, buildings, id_name, unique_id):
             single_geom.append(p)
 
     blocks_gdf = gpd.GeoDataFrame(geometry=gpd.GeoSeries(single_geom))
-    weights_matrix = Queen.from_dataframe(blocks_gdf, silence_warnings=True)
+    spatial_weights = Queen.from_dataframe(blocks_gdf, silence_warnings=True)
 
     patches = {}
     jID = 1
@@ -594,14 +594,14 @@ def blocks(cells, streets, buildings, id_name, unique_id):
         else:
             to_join = [idx]  # list of indices which should be joined together
             neighbours = []  # list of neighbours
-            weights = weights_matrix.neighbors[idx]  # neighbours from spatial weights
+            weights = spatial_weights.neighbors[idx]  # neighbours from spatial weights
             for w in weights:
                 neighbours.append(w)  # make a list from weigths
 
             for n in neighbours:
                 while n not in to_join:  # until there is some neighbour which is not in to_join
                     to_join.append(n)
-                    weights = weights_matrix.neighbors[n]
+                    weights = spatial_weights.neighbors[n]
                     for w in weights:
                         neighbours.append(w)  # extend neighbours by neighbours of neighbours :)
             for b in to_join:
diff --git a/momepy/intensity.py b/momepy/intensity.py
index ea49813..f48926c 100644
--- a/momepy/intensity.py
+++ b/momepy/intensity.py
@@ -9,96 +9,6 @@ import pandas as pd
 import collections
 
 
-def radius(gpd_df, cpt, radius):
-    """
-    Get a list of indices of objects within radius.
-
-    Parameters
-    ----------
-    gpd_df : GeoDataFrame
-        GeoDataFrame containing point objects to analyse
-    cpt : shapely.Point
-        shapely point representing the center of radius
-    radius : float
-        radius
-
-    Returns
-    -------
-    list
-        Return only the neighbour indices, sorted by distance in ascending order
-
-    Notes
-    ---------
-    https://stackoverflow.com/questions/44622233/rtree-count-points-in-the-neighbourhoods-within-each-point-of-another-set-of-po
-
-    """
-    # Spatial index
-    sindex = gpd_df.sindex
-    # Bounding box of rtree search (West, South, East, North)
-    bbox = (cpt.x - radius, cpt.y - radius, cpt.x + radius, cpt.y + radius)
-    # Potential neighbours
-    good = []
-    for n in sindex.intersection(bbox):
-        dist = cpt.distance(gpd_df['geometry'].iloc[n])
-        if dist < radius:
-            good.append((dist, n))
-    # Sort list in ascending order by `dist`, then `n`
-    good.sort()
-    # Return only the neighbour indices, sorted by distance in ascending order
-    return [x[1] for x in good]
-
-
-def frequency(objects, look_for, id_column='uID', rad=400):
-    """
-    Calculate frequency (count) of objects in a given radius.
-
-    .. math::
-        count
-
-    Parameters
-    ----------
-    objects : GeoDataFrame
-        GeoDataFrame containing objects to analyse
-    look_for : GeoDataFrame
-        GeoDataFrame with measured objects (could be the same as objects)
-    id_column : str
-        name of the column with unique id
-
-    Returns
-    -------
-    Series
-        Series containing resulting values.
-
-    References
-    ---------
-
-    Notes
-    -----
-    Might be faster using libpysal DistanceBand
-
-    """
-
-    print('Calculating frequency...')
-
-    objects_centroids = objects.copy()
-    objects_centroids['geometry'] = objects_centroids.centroid
-
-    look_for_centroids = look_for.copy()
-    look_for_centroids['geometry'] = look_for_centroids.centroid
-
-    # define empty list for results
-    results_list = []
-
-    for index, row in tqdm(objects_centroids.iterrows(), total=objects_centroids.shape[0]):
-        neighbours = radius(look_for_centroids, row['geometry'], rad)
-        results_list.append(len(neighbours))
-
-    series = pd.Series(results_list)
-
-    print('Frequency calculated.')
-    return series
-
-
 def covered_area_ratio(objects, look_for, area_column, look_for_area_column, id_column="uID"):
     """
     Calculate covered area ratio of objects.
@@ -250,7 +160,7 @@ def elements_in_block(blocks, elements, left_id, right_id, weighted=False):
     return series
 
 
-def courtyards(objects, block_id, weights_matrix=None):
+def courtyards(objects, block_id, spatial_weights=None):
     """
     Calculate the number of courtyards within the joined structure.
 
@@ -260,7 +170,7 @@ def courtyards(objects, block_id, weights_matrix=None):
         GeoDataFrame containing objects to analyse
     block_id : str
         name of the column where is stored block ID
-    weights_matrix : libpysal.weights, optional
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix will be calculated
         based on objects. It is to denote adjacent buildings.
 
@@ -282,10 +192,10 @@ def courtyards(objects, block_id, weights_matrix=None):
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
     # if weights matrix is not passed, generate it from objects
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Calculating spatial weights...')
         from libpysal.weights import Queen
-        weights_matrix = Queen.from_dataframe(objects, silence_warnings=True)
+        spatial_weights = Queen.from_dataframe(objects, silence_warnings=True)
         print('Spatial weights ready...')
 
     # dict to store nr of courtyards for each uID
@@ -298,14 +208,14 @@ def courtyards(objects, block_id, weights_matrix=None):
         else:
             to_join = [index]  # list of indices which should be joined together
             neighbours = []  # list of neighbours
-            weights = weights_matrix.neighbors[index]  # neighbours from spatial weights
+            weights = spatial_weights.neighbors[index]  # neighbours from spatial weights
             for w in weights:
                 neighbours.append(w)  # make a list from weigths
 
             for n in neighbours:
                 while n not in to_join:  # until there is some neighbour which is not in to_join
                     to_join.append(n)
-                    weights = weights_matrix.neighbors[n]
+                    weights = spatial_weights.neighbors[n]
                     for w in weights:
                         neighbours.append(w)  # extend neighbours by neighbours of neighbours :)
             joined = objects.iloc[to_join]
@@ -325,7 +235,7 @@ def courtyards(objects, block_id, weights_matrix=None):
     return series
 
 
-def gross_density(objects, buildings, area, character, weights_matrix=None, order=3, unique_id='uID'):
+def gross_density(objects, buildings, area, character, spatial_weights=None, order=3, unique_id='uID'):
     """
     Calculate the density
 
@@ -342,7 +252,7 @@ def gross_density(objects, buildings, area, character, weights_matrix=None, orde
         name of the column with area values
     character : str
         name of the column with values of target character for density calculation
-    weights_matrix : libpysal.weights, optional
+    spatial_weights : libpysal.weights, optional
         spatial weights matrix - If None, Queen contiguity matrix of selected order will be calculated
         based on objects.
     order : int
@@ -371,15 +281,15 @@ def gross_density(objects, buildings, area, character, weights_matrix=None, orde
     if not all(objects.index == range(len(objects))):
         raise ValueError('Index is not consecutive range 0:x, spatial weights will not match objects.')
 
-    if weights_matrix is None:
+    if spatial_weights is None:
         print('Generating weights matrix (Queen) of {} topological steps...'.format(order))
         from momepy import Queen_higher
         # matrix to define area of analysis (more steps)
-        weights_matrix = Queen_higher(k=order, geodataframe=objects)
+        spatial_weights = Queen_higher(k=order, geodataframe=objects)
 
     # iterating over rows one by one
     for index, row in tqdm(objects.iterrows(), total=objects.shape[0]):
-        neighbours_id = weights_matrix.neighbors[index]
+        neighbours_id = spatial_weights.neighbors[index]
         neighbours_id.append(index)
         neighbours = objects.iloc[neighbours_id]
 
diff --git a/momepy/utils.py b/momepy/utils.py
index 383ac86..3ec1dd0 100644
--- a/momepy/utils.py
+++ b/momepy/utils.py
@@ -184,17 +184,17 @@ def multi2single(gdf):
     return gdf_singlepoly
 
 
-def limit_range(vals, mode):
+def limit_range(vals, rng):
     """
-    Extract values within selected mode ('id', 'iq')
+    Extract values within selected range
 
     Parameters
     ----------
     vals : array
 
-    mode : str
-        'iq' for interquartile range, 'id' for interdecile range
-
+    rng : Two-element sequence containing floats in range of [0,100], optional
+        Percentiles over which to compute the range. Each must be
+        between 0 and 100, inclusive. The order of the elements is not important.
     Returns
     -------
     array
@@ -202,14 +202,9 @@ def limit_range(vals, mode):
     """
     if len(vals) > 2:
         limited = []
-        if mode == 'iq':
-            mval = 25
-            xval = 75
-        elif mode == 'id':
-            mval = 10
-            xval = 90
-        lower = np.percentile(vals, mval)
-        higher = np.percentile(vals, xval)
+        rng = sorted(rng)
+        lower = np.percentile(vals, rng[0])
+        higher = np.percentile(vals, rng[1])
         for x in vals:
             if x >= lower and x <= higher:
                 limited.append(x)


Test patch:
diff --git a/tests/test_dimension.py b/tests/test_dimension.py
index 90f81a3..47678a6 100644
--- a/tests/test_dimension.py
+++ b/tests/test_dimension.py
@@ -105,8 +105,8 @@ class TestDimensions:
         self.df_tessellation['area'] = area = self.df_tessellation.geometry.area
         self.df_tessellation['mesh_ar'] = mm.mean_character(self.df_tessellation, spatial_weights, values='area')
         self.df_tessellation['mesh_array'] = mm.mean_character(self.df_tessellation, spatial_weights, values=area)
-        self.df_tessellation['mesh_id'] = mm.mean_character(self.df_tessellation, spatial_weights, values='area', mode='id')
-        self.df_tessellation['mesh_iq'] = mm.mean_character(self.df_tessellation, spatial_weights, values='area', mode='iq')
+        self.df_tessellation['mesh_id'] = mm.mean_character(self.df_tessellation, spatial_weights, values='area', rng=(10, 90))
+        self.df_tessellation['mesh_iq'] = mm.mean_character(self.df_tessellation, spatial_weights, values='area', rng=(25, 75))
         neighbours = spatial_weights.neighbors[38]
         total_area = sum(self.df_tessellation.iloc[neighbours].geometry.area) + self.df_tessellation.geometry.area[38]
         check = total_area / (len(neighbours) + 1)
@@ -158,3 +158,17 @@ class TestDimensions:
         sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
         weighted = mm.weighted_character(self.df_buildings, self.df_tessellation, 'height', 'uID', sw, area)
         assert weighted[38] == 18.301521351817303
+
+    def test_covered_area(self):
+        sw = Queen_higher(geodataframe=self.df_tessellation, k=1)
+        covered = mm.covered_area(self.df_tessellation)
+        covered_sw = mm.covered_area(self.df_tessellation, sw)
+        assert covered[0] == covered_sw[0]
+        assert covered[0] == 24115.66721833942
+
+    def test_wall(self):
+        sw = Queen_higher(geodataframe=self.df_buildings, k=1)
+        wall = mm.wall(self.df_buildings)
+        wall_sw = mm.wall(self.df_buildings, sw)
+        assert wall[0] == wall_sw[0]
+        assert wall[0] == 137.2106961418436
diff --git a/tests/test_distribution.py b/tests/test_distribution.py
index e5589a3..92ed0dc 100644
--- a/tests/test_distribution.py
+++ b/tests/test_distribution.py
@@ -42,15 +42,15 @@ class TestDistribution:
         assert self.df_buildings['street_a_arr'][0] == check
 
     def test_cell_alignment(self):
-        self.df_buildings['orient'] = orient = mm.orientation(self.df_buildings)
-        self.df_tessellation['orient'] = tess_orient = mm.orientation(self.df_tessellation)
+        self.df_buildings['orient'] = mm.orientation(self.df_buildings)
+        self.df_tessellation['orient'] = mm.orientation(self.df_tessellation)
         self.df_buildings['c_align'] = mm.cell_alignment(self.df_buildings, self.df_tessellation, 'orient', 'orient', 'uID')
         check = abs(self.df_buildings['orient'][0] -
                     self.df_tessellation[self.df_tessellation['uID'] == self.df_buildings['uID'][0]]['orient'].iloc[0])
         assert self.df_buildings['c_align'][0] == check
 
     def test_alignment(self):
-        self.df_buildings['orient'] = orient = mm.orientation(self.df_buildings)
+        self.df_buildings['orient'] = mm.orientation(self.df_buildings)
         self.df_buildings['align'] = mm.alignment(self.df_buildings, 'orient', self.df_tessellation, 'uID')
         sw = Queen.from_dataframe(self.df_tessellation)
         self.df_buildings['align_sw'] = mm.alignment(self.df_buildings, 'orient', self.df_tessellation, 'uID', sw)
@@ -84,7 +84,7 @@ class TestDistribution:
         self.df_buildings['adj'] = mm.building_adjacency(self.df_buildings, self.df_tessellation)
         sw = Queen.from_dataframe(self.df_buildings)
         swh = mm.Queen_higher(k=3, geodataframe=self.df_tessellation)
-        self.df_buildings['adj_sw'] = mm.building_adjacency(self.df_buildings, self.df_tessellation, weights_matrix=sw, weights_matrix_higher=swh)
+        self.df_buildings['adj_sw'] = mm.building_adjacency(self.df_buildings, self.df_tessellation, spatial_weights=sw, spatial_weights_higher=swh)
         check = 0.2613824113909074
         assert self.df_buildings['adj'].mean() == check
         assert self.df_buildings['adj_sw'].mean() == check
diff --git a/tests/test_diversity.py b/tests/test_diversity.py
new file mode 100644
index 0000000..c2b0e97
--- /dev/null
+++ b/tests/test_diversity.py
@@ -0,0 +1,48 @@
+import pytest
+import momepy as mm
+import geopandas as gpd
+import numpy as np
+
+from momepy import Queen_higher
+
+
+class TestDiversity:
+
+    def setup_method(self):
+
+        test_file_path = mm.datasets.get_path('bubenec')
+        self.df_buildings = gpd.read_file(test_file_path, layer='buildings')
+        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., 30., 144)
+        self.df_tessellation['area'] = mm.area(self.df_tessellation)
+
+    def test_rng(self):
+        full = mm.rng(self.df_tessellation, 'area', order=3)
+        assert full[0] == 8255.372874447059
+        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
+        full_sw = mm.rng(self.df_tessellation, 'area', spatial_weights=sw)
+        assert full_sw[0] == 8255.372874447059
+        area = self.df_tessellation['area']
+        full2 = mm.rng(self.df_tessellation, area, order=3)
+        assert full2[0] == 8255.372874447059
+        limit = mm.rng(self.df_tessellation, 'area', spatial_weights=sw, rng=(10, 90))
+        assert limit[0] == 4122.139212736442
+
+    def test_theil(self):
+        full = mm.theil(self.df_tessellation, 'area', order=3)
+        assert full[0] == 0.2574468431886534
+        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
+        full_sw = mm.theil(self.df_tessellation, 'area', spatial_weights=sw)
+        assert full_sw[0] == 0.2574468431886534
+        limit = mm.theil(self.df_tessellation, 'area', spatial_weights=sw, rng=(10, 90))
+        assert limit[0] == 0.13302952097969373
+
+    def test_simpson(self):
+        ht = mm.simpson(self.df_tessellation, 'area', order=3)
+        assert ht[0] == 0.385
+        sw = Queen_higher(k=3, geodataframe=self.df_tessellation)
+        ht_sw = mm.simpson(self.df_tessellation, 'area', spatial_weights=sw)
+        assert ht_sw[0] == 0.385
+        quan_sw = mm.simpson(self.df_tessellation, 'area', spatial_weights=sw, binning='quantiles', k=3)
+        assert quan_sw[0] == 0.395
diff --git a/tests/test_intensity.py b/tests/test_intensity.py
index 010797a..bba6017 100644
--- a/tests/test_intensity.py
+++ b/tests/test_intensity.py
@@ -47,7 +47,7 @@ class TestDistribution:
     def test_gross_density(self):
         dens = mm.gross_density(self.df_tessellation, self.df_buildings, 'area', 'fl_area')
         sw = mm.Queen_higher(k=3, geodataframe=self.df_tessellation)
-        dens2 = mm.gross_density(self.df_tessellation, self.df_buildings, 'area', 'fl_area', weights_matrix=sw)
+        dens2 = mm.gross_density(self.df_tessellation, self.df_buildings, 'area', 'fl_area', spatial_weights=sw)
         check = 1.6615871155383324
         assert dens.mean() == check
         assert dens2.mean() == check
diff --git a/tests/test_utils.py b/tests/test_utils.py
index ba5a39f..d00e3ff 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -69,8 +69,8 @@ class TestUtils:
         assert len(nodes) == 32
 
     def test_limit_range(self):
-        assert mm.limit_range(range(10), mode='iq') == [3, 4, 5, 6]
-        assert mm.limit_range(range(10), mode='id') == [1, 2, 3, 4, 5, 6, 7, 8]
+        assert mm.limit_range(range(10), rng=(25, 75)) == [3, 4, 5, 6]
+        assert mm.limit_range(range(10), rng=(10, 90)) == [1, 2, 3, 4, 5, 6, 7, 8]
 
     def test_preprocess(self):
         test_file_path2 = mm.datasets.get_path('tests')