text: #!/usr/bin/env python
  # -*- coding: utf-8 -*-
  
  # shape.py
  # definitions of shape characters
  
  import math
  import random
  
  import numpy as np
  import pandas as pd
  from shapely.geometry import Point
  from tqdm import tqdm  # progress bar
  
  __all__ = [
      "FormFactor",
      "FractalDimension",
      "VolumeFacadeRatio",
      "CircularCompactness",
      "SquareCompactness",
      "Convexity",
      "CourtyardIndex",
      "Rectangularity",
      "ShapeIndex",
      "Corners",
      "Squareness",
      "EquivalentRectangularIndex",
      "Elongation",
      "CentroidCorners",
      "Linearity",
      "CompactnessWeightedAxis",
  ]
  
  
  class FormFactor:
      """
      Calculates form factor of each object in given GeoDataFrame.
  
      .. math::
          area \\over {volume^{2 \\over 3}}
  
      Adapted from :cite:`bourdic2012`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      volumes : str, list, np.array, pd.Series
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored volume value.
          (To calculate volume you can use :py:func:`momepy.volume`)
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      volumes : Series
          Series containing used volume values
      areas : Series
          Series containing used area values
  
      Examples
      --------
      >>> buildings_df['formfactor'] = momepy.FormFactor(buildings_df, 'volume').series
      >>> buildings_df.formfactor[0]
      1.9385988170288635
  
      >>> volume = momepy.Volume(buildings_df, 'height').series
      >>> buildings_df['formfactor'] = momepy.FormFactor(buildings_df, volume).series
      >>> buildings_df.formfactor[0]
      1.9385988170288635
  
      """
  
      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"
          self.areas = gdf[areas]
          zeros = gdf[volumes] == 0
          res = np.empty(len(gdf))
          res[zeros] = 0
          res[~zeros] = gdf[areas][~zeros] / (gdf[volumes][~zeros] ** (2 / 3))
          self.series = pd.Series(res, index=gdf.index)
  
  
  class FractalDimension:
      """
      Calculates fractal dimension of each object in given GeoDataFrame.
  
      .. math::
          {2log({{perimeter} \\over {4}})} \\over log(area)
  
      Based on :cite:`mcgarigal1995fragstats`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
      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. If set to ``None``, function will calculate perimeters
          during the process without saving them separately.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      perimeters : Series
          Series containing used perimeter values
      areas : Series
          Series containing used area values
  
      Examples
      --------
      >>> buildings_df['fractal'] = momepy.FractalDimension(buildings_df,
      ...                                                   'area',
      ...                                                   'peri').series
      >>> buildings_df.fractal[0]
      1.0726778567038908
      """
  
      def __init__(self, gdf, areas=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 areas is None:
              areas = gdf.geometry.area
          if not isinstance(areas, str):
              gdf["mm_a"] = areas
              areas = "mm_a"
  
          self.series = pd.Series(
              (2 * np.log(gdf[perimeters] / 4)) / np.log(gdf[areas]), index=gdf.index
          )
  
  
  class VolumeFacadeRatio:
      """
      Calculates volume/facade ratio of each object in given GeoDataFrame.
  
      .. math::
          volume \\over perimeter * height
  
      Adapted from :cite:`schirmer2015`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      heights : str, list, np.array, pd.Series (default None)
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored height value
      volumes : str, list, np.array, pd.Series (default None)
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored volume value
      perimeters : , list, np.array, pd.Series (default None)
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored perimeter value
  
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      perimeters : Series
          Series containing used perimeter values
      volumes : Series
          Series containing used volume values
  
      Examples
      --------
      >>> buildings_df['vfr'] = momepy.VolumeFacadeRatio(buildings_df, 'height').series
      >>> buildings_df.vfr[0]
      5.310715735236504
      """
  
      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
                  volumes = "mm_v"
          self.volumes = gdf[volumes]
  
          self.series = gdf[volumes] / (gdf[perimeters] * gdf[heights])
  
  
  # Smallest enclosing circle - Library (Python)
  
  # Copyright (c) 2017 Project Nayuki
  # https://www.nayuki.io/page/smallest-enclosing-circle
  
  # This program is free software: you can redistribute it and/or modify
  # it under the terms of the GNU Lesser General Public License as published by
  # the Free Software Foundation, either version 3 of the License, or
  # (at your option) any later version.
  
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU Lesser General Public License for more details.
  
  # You should have received a copy of the GNU Lesser General Public License
  # along with this program (see COPYING.txt and COPYING.LESSER.txt).
  # If not, see <http://www.gnu.org/licenses/>.
  
  # Data conventions: A point is a pair of floats (x, y).
  # A circle is a triple of floats (center x, center y, radius).
  
  # Returns the smallest circle that encloses all the given points.
  # Runs in expected O(n) time, randomized.
  # Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
  
  # Output: A triple of floats representing a circle.
  # Note: If 0 points are given, None is returned. If 1 point is given,
  # a circle of radius 0 is returned.
  #
  # Initially: No boundary points known
  
  
  def _make_circle(points):
      # Convert to float and randomize order
      shuffled = [(float(x), float(y)) for (x, y) in points]
      random.shuffle(shuffled)
  
      # Progressively add points to circle or recompute circle
      c = None
      for (i, p) in enumerate(shuffled):
          if c is None or not _is_in_circle(c, p):
              c = _make_circle_one_point(shuffled[: i + 1], p)
      return c
  
  
  # One boundary point known
  def _make_circle_one_point(points, p):
      c = (p[0], p[1], 0.0)
      for (i, q) in enumerate(points):
          if not _is_in_circle(c, q):
              if c[2] == 0.0:
                  c = _make_diameter(p, q)
              else:
                  c = _make_circle_two_points(points[: i + 1], p, q)
      return c
  
  
  # Two boundary points known
  def _make_circle_two_points(points, p, q):
      circ = _make_diameter(p, q)
      left = None
      right = None
      px, py = p
      qx, qy = q
  
      # For each point not in the two-point circle
      for r in points:
          if _is_in_circle(circ, r):
              continue
  
          # Form a circumcircle and classify it on left or right side
          cross = _cross_product(px, py, qx, qy, r[0], r[1])
          c = _make_circumcircle(p, q, r)
          if c is None:
              continue
          elif cross > 0.0 and (
              left is None
              or _cross_product(px, py, qx, qy, c[0], c[1])
              > _cross_product(px, py, qx, qy, left[0], left[1])
          ):
              left = c
          elif cross < 0.0 and (
              right is None
              or _cross_product(px, py, qx, qy, c[0], c[1])
              < _cross_product(px, py, qx, qy, right[0], right[1])
          ):
              right = c
  
      # Select which circle to return
      if left is None and right is None:
          return circ
      if left is None:
          return right
      if right is None:
          return left
      if left[2] <= right[2]:
          return left
      return right
  
  
  def _make_circumcircle(p0, p1, p2):
      # Mathematical algorithm from Wikipedia: Circumscribed circle
      ax, ay = p0
      bx, by = p1
      cx, cy = p2
      ox = (min(ax, bx, cx) + max(ax, bx, cx)) / 2.0
      oy = (min(ay, by, cy) + max(ay, by, cy)) / 2.0
      ax -= ox
      ay -= oy
      bx -= ox
      by -= oy
      cx -= ox
      cy -= oy
      d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0
      if d == 0.0:
          return None
      x = (
          ox
          + (
              (ax * ax + ay * ay) * (by - cy)
              + (bx * bx + by * by) * (cy - ay)
              + (cx * cx + cy * cy) * (ay - by)
          )
          / d
      )
      y = (
          oy
          + (
              (ax * ax + ay * ay) * (cx - bx)
              + (bx * bx + by * by) * (ax - cx)
              + (cx * cx + cy * cy) * (bx - ax)
          )
          / d
      )
      ra = math.hypot(x - p0[0], y - p0[1])
      rb = math.hypot(x - p1[0], y - p1[1])
      rc = math.hypot(x - p2[0], y - p2[1])
      return (x, y, max(ra, rb, rc))
  
  
  def _make_diameter(p0, p1):
      cx = (p0[0] + p1[0]) / 2.0
      cy = (p0[1] + p1[1]) / 2.0
      r0 = math.hypot(cx - p0[0], cy - p0[1])
      r1 = math.hypot(cx - p1[0], cy - p1[1])
      return (cx, cy, max(r0, r1))
  
  
  _MULTIPLICATIVE_EPSILON = 1 + 1e-14
  
  
  def _is_in_circle(c, p):
      return (
          c is not None
          and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON
      )
  
  
  # Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2).
  def _cross_product(x0, y0, x1, y1, x2, y2):
      return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
  
  
  # end of Nayuiki script to define the smallest enclosing circle
  
  
  # calculate the area of circumcircle
  def _circle_area(points):
      if len(points[0]) == 3:
          points = [x[:2] for x in points]
      circ = _make_circle(points)
      return math.pi * circ[2] ** 2
  
  
  def _circle_radius(points):
      if len(points[0]) == 3:
          points = [x[:2] for x in points]
      circ = _make_circle(points)
      return circ[2]
  
  
  class CircularCompactness:
      """
      Calculates compactness index of each object in given GeoDataFrame.
  
      .. math::
          area \\over \\textit{area of enclosing circle}
  
      Adapted from :cite:`dibble2017`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      areas : Series
          Series containing used area values
  
      Examples
      --------
      >>> buildings_df['comp'] = momepy.CircularCompactness(buildings_df, 'area').series
      >>> buildings_df['comp'][0]
      0.572145421828038
      """
  
      def __init__(self, gdf, areas=None):
          self.gdf = gdf
  
          if areas is None:
              areas = gdf.geometry.area
          elif isinstance(areas, str):
              areas = gdf[areas]
          self.areas = areas
          hull = gdf.convex_hull.exterior
          radius = hull.apply(
              lambda g: _circle_radius(list(g.coords)) if g is not None else None
          )
          self.series = areas / (np.pi * radius ** 2)
  
  
  class SquareCompactness:
      """
      Calculates compactness index of each object in given GeoDataFrame.
  
      .. math::
          \\begin{equation*}
          \\left(\\frac{4 \\sqrt{area}}{perimeter}\\right) ^ 2
          \\end{equation*}
  
      Adapted from :cite:`feliciotti2018`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
      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. If set to ``None``, function will calculate perimeters
          during the process without saving them separately.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      areas : Series
          Series containing used area values
      perimeters : Series
          Series containing used perimeter values
  
      Examples
      --------
      >>> buildings_df['squ_comp'] = momepy.SquareCompactness(buildings_df).series
      >>> buildings_df['squ_comp'][0]
      0.6193872538650996
  
      """
  
      def __init__(self, gdf, areas=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 areas is None:
              areas = gdf.geometry.area
          if not isinstance(areas, str):
              gdf["mm_a"] = areas
              areas = "mm_a"
          self.areas = gdf[areas]
          self.series = ((np.sqrt(gdf[areas]) * 4) / gdf[perimeters]) ** 2
  
  
  class Convexity:
      """
      Calculates Convexity index of each object in given GeoDataFrame.
  
      .. math::
          area \\over \\textit{convex hull area}
  
      Adapted from :cite:`dibble2017`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      areas : Series
          Series containing used area values
  
      Examples
      --------
      >>> buildings_df['convexity'] = momepy.Convexity(buildings_df).series
      >>> buildings_df['convexity'][0]
      0.8151964258521672
      """
  
      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"
          self.areas = gdf[areas]
          self.series = gdf[areas] / gdf.geometry.convex_hull.area
  
  
  class CourtyardIndex:
      """
      Calculates courtyard index of each object in given GeoDataFrame.
  
      .. math::
          \\textit{area of courtyards} \\over \\textit{total area}
  
      Adapted from :cite:`schirmer2015`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      courtyard_areas : str, list, np.array, pd.Series
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored area value
          (To calculate volume you can use :py:class:`momepy.CourtyardArea`)
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : 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
  
      Examples
      --------
      >>> buildings_df['courtyard_index'] = momepy.CourtyardIndex(buildings,
      ...                                                         'courtyard_area',
      ...                                                         'area').series
      >>> buildings_df.courtyard_index[80]
      0.16605915738643523
      """
  
      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"
          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"
          self.areas = gdf[areas]
          self.series = gdf[courtyard_areas] / gdf[areas]
  
  
  class Rectangularity:
      """
      Calculates rectangularity of each object in given GeoDataFrame.
  
      .. math::
          {area \\over \\textit{minimum bounding rotated rectangle area}}
  
      Adapted from :cite:`dibble2017`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      areas : Series
          Series containing used area values
  
      Examples
      --------
      >>> buildings_df['rect'] = momepy.Rectangularity(buildings_df, 'area').series
      100%|██████████| 144/144 [00:00<00:00, 866.62it/s]
      >>> buildings_df.rect[0]
      0.6942676157646379
      """
  
      def __init__(self, gdf, areas=None):
          # TODO: vectorize minimum_rotated_rectangle after pygeos implementation
          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"
          self.areas = gdf[areas]
          self.series = gdf.apply(
              lambda row: row[areas] / (row.geometry.minimum_rotated_rectangle.area),
              axis=1,
          )
  
  
  class ShapeIndex:
      """
      Calculates shape index of each object in given GeoDataFrame.
  
      .. math::
          {\\sqrt{{area} \\over {\\pi}}} \\over {0.5 * \\textit{longest axis}}
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      longest_axis : str, list, np.array, pd.Series
          the name of the dataframe column, ``np.array``, or ``pd.Series`` where is
          stored longest axis value
      areas : str, list, np.array, pd.Series (default 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.
  
      Attributes
      ----------
      series : 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
  
      Examples
      --------
      >>> buildings_df['shape_index'] = momepy.ShapeIndex(buildings_df,
      ...                                                 longest_axis='long_ax',
      ...                                                 areas='area').series
      100%|██████████| 144/144 [00:00<00:00, 5558.33it/s]
      >>> buildings_df['shape_index'][0]
      0.7564029493781987
      """
  
      def __init__(self, gdf, longest_axis, areas=None):
          self.gdf = gdf
          gdf = gdf.copy()
  
          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"
          self.areas = gdf[areas]
          self.series = pd.Series(
              np.sqrt(gdf[areas] / np.pi) / (0.5 * gdf[longest_axis]), index=gdf.index
          )
  
  
  class Corners:
      """
      Calculates number of corners of each object in given GeoDataFrame.
  
      Uses only external shape (``shapely.geometry.exterior``), courtyards are not
      included.
  
      .. math::
          \\sum corner
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      verbose : bool (default True)
          if True, shows progress bars in loops and indication of steps
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
  
  
      Examples
      --------
      >>> buildings_df['corners'] = momepy.Corners(buildings_df).series
      100%|██████████| 144/144 [00:00<00:00, 1042.15it/s]
      >>> buildings_df.corners[0]
      24
  
  
      """
  
      def __init__(self, gdf, verbose=True):
          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)
  
              # TODO: add arg to specify these values
              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 geom in tqdm(gdf.geometry, total=gdf.shape[0], disable=not verbose):
              if geom.type == "Polygon":
                  corners = 0  # define empty variables
                  points = list(geom.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:
                              continue
  
                      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
              elif geom.type == "MultiPolygon":
                  corners = 0  # define empty variables
                  for g in geom.geoms:
                      points = list(g.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:
                                  continue
  
                          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
              else:
                  corners = np.nan
  
              results_list.append(corners)
  
          self.series = pd.Series(results_list, index=gdf.index)
  
  
  class Squareness:
      """
      Calculates squareness of each object in given GeoDataFrame.
  
      Uses only external shape (``shapely.geometry.exterior``), courtyards are not
      included.
  
      .. math::
          \\mu=\\frac{\\sum_{i=1}^{N} d_{i}}{N}
  
      where :math:`d` is the deviation of angle of corner :math:`i` from 90 degrees.
  
      Adapted from :cite:`dibble2017`.
  
      Returns ``np.nan`` for MultiPolygons.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      verbose : bool (default True)
          if True, shows progress bars in loops and indication of steps
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
  
      Examples
      --------
      >>> buildings_df['squareness'] = momepy.Squareness(buildings_df).series
      100%|██████████| 144/144 [00:01<00:00, 129.49it/s]
      >>> buildings_df.squareness[0]
      3.7075816043359864
      """
  
      def __init__(self, gdf, verbose=True):
          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 geom in tqdm(gdf.geometry, total=gdf.shape[0], disable=not verbose):
              if geom.type == "Polygon":
                  angles = []
                  points = list(geom.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:
                              continue
                  deviations = [abs(90 - i) for i in angles]
                  results_list.append(np.mean(deviations))
  
              else:
                  results_list.append(np.nan)
  
          self.series = pd.Series(results_list, index=gdf.index)
  
  
  class EquivalentRectangularIndex:
      """
      Calculates equivalent rectangular index of each object in given GeoDataFrame.
  
      .. math::
          \\sqrt{{area} \\over \\textit{area of bounding rectangle}} *
          {\\textit{perimeter of bounding rectangle} \\over {perimeter}}
  
      Based on :cite:`basaraner2017`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
      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. If set to ``None``, function will calculate perimeters
          during the process without saving them separately.
  
      Attributes
      ----------
      series : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
      areas : Series
          Series containing used area values
      perimeters : Series
          Series containing used perimeter values
  
      Examples
      --------
      >>> buildings_df['eri'] = momepy.EquivalentRectangularIndex(buildings_df,
      ...                                                         'area',
      ...                                                         'peri').series
      >>> buildings_df['eri'][0]
      0.7879229963118455
      """
  
      def __init__(self, gdf, areas=None, perimeters=None):
          self.gdf = gdf
          # define empty list for results
  
          if perimeters is None:
              perimeters = gdf.geometry.length
          else:
              if isinstance(perimeters, str):
                  perimeters = gdf[perimeters]
  
          self.perimeters = perimeters
  
          if areas is None:
              areas = gdf.geometry.area
          else:
              if isinstance(areas, str):
                  areas = gdf[areas]
  
          self.areas = areas
          # TODO: vectorize minimum_rotated_rectangle after pygeos implementation
          bbox = gdf.geometry.apply(lambda g: g.minimum_rotated_rectangle)
          res = np.sqrt(areas / bbox.area) * (bbox.length / perimeters)
  
          self.series = pd.Series(res, index=gdf.index)
  
  
  class Elongation:
      """
      Calculates elongation of object seen as elongation of
      its minimum bounding rectangle.
  
      .. math::
          {{p - \\sqrt{p^2 - 16a}} \\over {4}} \\over
          {{{p} \\over {2}} - {{p - \\sqrt{p^2 - 16a}} \\over {4}}}
  
      where `a` is the area of the object and `p` its perimeter.
  
      Based on :cite:`gil2012`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
  
      Attributes
      ----------
      e : Series
          Series containing resulting values
      gdf : GeoDataFrame
          original GeoDataFrame
  
      Examples
      --------
      >>> buildings_df['elongation'] = momepy.Elongation(buildings_df).series
      >>> buildings_df['elongation'][0]
      0.9082437463675544
      """
  
      def __init__(self, gdf):
          self.gdf = gdf
  
          # TODO: vectorize minimum_rotated_rectangle after pygeos implementation
          bbox = gdf.geometry.apply(lambda g: g.minimum_rotated_rectangle)
          a = bbox.area
          p = bbox.length
          cond1 = p ** 2
          cond2 = 16 * a
          bigger = cond1 >= cond2
          sqrt = np.empty(len(a))
          sqrt[bigger] = cond1[bigger] - cond2[bigger]
          sqrt[~bigger] = 0
  
          # calculate both width/length and length/width
          elo1 = ((p - np.sqrt(sqrt)) / 4) / ((p / 2) - ((p - np.sqrt(sqrt)) / 4))
          elo2 = ((p + np.sqrt(sqrt)) / 4) / ((p / 2) - ((p + np.sqrt(sqrt)) / 4))
  
          # use the smaller one (e.g. shorter/longer)
          res = np.empty(len(a))
          res[elo1 <= elo2] = elo1[elo1 <= elo2]
          res[~(elo1 <= elo2)] = elo2[~(elo1 <= elo2)]
  
          self.series = pd.Series(res, index=gdf.index)
  
  
  class CentroidCorners:
      """
      Calculates mean distance centroid - corners and st. deviation.
  
      .. math::
          \\overline{x}=\\frac{1}{n}\\left(\\sum_{i=1}^{n} dist_{i}\\right);
          \\space \\mathrm{SD}=\\sqrt{\\frac{\\sum|x-\\overline{x}|^{2}}{n}}
  
      Adapted from :cite:`schirmer2015` and :cite:`cimburova2017`.
  
      Returns ``np.nan`` for MultiPolygons.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      verbose : bool (default True)
          if True, shows progress bars in loops and indication of steps
  
      Attributes
      ----------
      mean : Series
          Series containing mean distance values.
      std : Series
          Series containing standard deviation values.
      gdf : GeoDataFrame
          original GeoDataFrame
  
      Examples
      --------
      >>> ccd = momepy.CentroidCorners(buildings_df)
      100%|██████████| 144/144 [00:00<00:00, 846.58it/s]
      >>> 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
      """
  
      def __init__(self, gdf, verbose=True):
          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 geom in tqdm(gdf.geometry, total=gdf.shape[0], disable=not verbose):
              if geom.type == "Polygon":
                  distances = []  # set empty list of distances
                  centroid = geom.centroid  # define centroid
                  points = list(geom.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:
                              continue
                  if not distances:  # circular buildings
                      if geom.has_z:
                          coords = [
                              (coo[0], coo[1]) for coo in geom.convex_hull.exterior.coords
                          ]
                      else:
                          coords = geom.convex_hull.exterior.coords
                      results_list.append(_circle_radius(coords))
                      results_list_sd.append(0)
                  else:
                      results_list.append(np.mean(distances))  # calculate mean
                      results_list_sd.append(np.std(distances))  # calculate st.dev
              else:
                  results_list.append(np.nan)
                  results_list_sd.append(np.nan)
  
          self.mean = pd.Series(results_list, index=gdf.index)
          self.std = pd.Series(results_list_sd, index=gdf.index)
  
  
  class Linearity:
      """
      Calculates linearity of each LineString object in given GeoDataFrame.
  
      .. math::
          \\frac{l_{euclidean}}{l_{segment}}
  
      where `l` is the length of the LineString. MultiLineString returns ``np.nan``.
  
      Adapted from :cite:`araldi2019`.
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      verbose : bool (default True)
          if True, shows progress bars in loops and indication of steps
  
      Attributes
      ----------
      series : Series
          Series containing mean distance values.
      gdf : GeoDataFrame
          original GeoDataFrame
  
      Examples
      --------
      >>> streets_df['linearity'] = momepy.Linearity(streets_df).series
      >>> streets_df['linearity'][0]
      1.0
      """
  
      def __init__(self, gdf, verbose=True):
          self.gdf = gdf
  
          euclidean = gdf.geometry.apply(
              lambda geom: self._dist(geom.coords[0], geom.coords[-1])
              if geom.type == "LineString"
              else np.nan
          )
          self.series = euclidean / gdf.geometry.length
  
      def _dist(self, a, b):
          return math.hypot(b[0] - a[0], b[1] - a[1])
  
  
  class CompactnessWeightedAxis:
      """
      Calculates compactness-weighted axis of each object in given GeoDataFrame.
  
      Initially designed for blocks.
  
      .. math::
          d_{i} \\times\\left(\\frac{4}{\\pi}-\\frac{16 (area_{i})}
          {perimeter_{i}^{2}}\\right)
  
      Parameters
      ----------
      gdf : GeoDataFrame
          GeoDataFrame containing objects
      areas : str, list, np.array, pd.Series (default 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.
      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. If set to ``None``, function will calculate perimeters
          during the process without saving them separately.
      longest_axis : str, list, np.array, pd.Series (default None)
          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.
  
      Attributes
      ----------
      series : 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.CompactnessWeightedAxis(blocks_df).series
      """
  
      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
                  perimeters = "mm_p"
          self.perimeters = gdf[perimeters]
          if longest_axis is None:
              from .dimension import LongestAxisLength
  
              gdf["mm_la"] = LongestAxisLength(gdf).series
              longest_axis = "mm_la"
          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"
          self.areas = gdf[areas]
          self.series = pd.Series(
              gdf[longest_axis]
              * ((4 / np.pi) - (16 * gdf[areas]) / ((gdf[perimeters]) ** 2)),
              index=gdf.index,
          )
  