docs: def area_of_polygon(polygon):
      """
      Returns the area of an OpenQuake polygon in square kilometres
      """
      lon0 = np.mean(polygon.lons)
      lat0 = np.mean(polygon.lats)
      # Transform to lamber equal area projection
      x, y = lonlat_to_laea(polygon.lons, polygon.lats, lon0, lat0)
      # Build shapely polygons
      poly = geometry.Polygon(zip(x, y))
      return poly.area