estimate-area command (#104)

* added estimate-area command

* added input tests

* tested out some feature input scenarios

* added 1cm flag

* edited some error messages for more clarity

* update Changelog and init for version 1.5.0

* moved over some tests and removed an exception caught

* fixed formatting for precision-testing.ldgeosjon

* updated min version for mercantile + supermercado

* edit 1cm error msg

CHANGELOG.md

@@ -1,5 +1,8 @@
 # Unreleased
 
+# 1.5.0 (2020-10-16)
+- Create estimate-area command 
+
 # 1.4.3 (2020-10-08)
 - Update Click version to 7.1.2 and fix assertions to pass all tests 
 

mapbox_tilesets/__init__.py

@@ -1,3 +1,3 @@
 """mapbox_tilesets package"""
 
-__version__ = "1.4.3"
+__version__ = "1.5.0"

mapbox_tilesets/scripts/cli.py

@@ -9,6 +9,8 @@
 
 import mapbox_tilesets
 from mapbox_tilesets import utils, errors
+from supermercado.super_utils import filter_features
+import builtins
 
 
 @click.version_option(version=mapbox_tilesets.__version__, message="%(version)s")
@@ -703,3 +705,70 @@ def list_sources(username, token=None):
             click.echo(source["id"])
     else:
         raise errors.TilesetsError(r.text)
+
+
+@cli.command("estimate-area")
+@cligj.features_in_arg
+@click.option(
+    "--precision",
+    "-p",
+    required=True,
+    type=click.Choice(["10m", "1m", "30cm", "1cm"]),
+    help="Precision level",
+)
+@click.option(
+    "--no-validation",
+    required=False,
+    is_flag=True,
+    help="Bypass source file validation",
+)
+@click.option(
+    "--force-1cm",
+    required=False,
+    is_flag=True,
+    help="Enables 1cm precision",
+)
+def estimate_area(features, precision, no_validation=False, force_1cm=False):
+    """Estimate area of features with a precision level.
+
+    tilesets estimate-area <features> <precision>
+
+    features must be a list of paths to local files containing GeoJSON feature collections or feature sequences from argument or stdin, or a list of string-encoded coordinate pairs of the form "[lng, lat]", or "lng, lat", or "lng lat".
+    """
+    area = 0
+    if precision == "1cm" and not force_1cm:
+        raise errors.TilesetsError(
+            "The --force-1cm flag must be present to enable 1cm precision area calculation and may take longer for large feature inputs or data with global extents. 1cm precision for tileset processing is only available upon request after contacting Mapbox support."
+        )
+    if precision != "1cm" and force_1cm:
+        raise errors.TilesetsError(
+            "The --force-1cm flag is enabled but the precision is not 1cm."
+        )
+
+    # builtins.list because there is a list command in the cli & will thrown an error
+    try:
+        features = builtins.list(filter_features(features))
+    except (ValueError, json.decoder.JSONDecodeError):
+        raise errors.TilesetsError(
+            "Error with feature parsing. Ensure that feature inputs are valid and formatted correctly. Try 'tilesets estimate-area --help' for help."
+        )
+    except Exception:
+        raise errors.TilesetsError("Error with feature filtering.")
+
+    # expect users to bypass source validation when users rerun command and their features passed validation previously
+    if not no_validation:
+        for feature in features:
+            utils.validate_geojson(feature)
+
+    area = utils.calculate_tiles_area(features, precision)
+    area = str(round(area))
+
+    click.echo(
+        json.dumps(
+            {
+                "km2": area,
+                "precision": precision,
+                "pricing_docs": "For more information, visit https://www.mapbox.com/pricing/#tilesets",
+            }
+        )
+    )

mapbox_tilesets/utils.py

@@ -1,8 +1,11 @@
 import os
 import re
 
+import numpy as np
+
 from jsonschema import validate
 from requests import Session
+from supermercado.burntiles import burn
 
 import mapbox_tilesets
 
@@ -103,3 +106,106 @@ def validate_geojson(feature):
     }
 
     return validate(instance=feature, schema=schema)
+
+
+def _convert_precision_to_zoom(precision):
+    """Converts precision to zoom level based on the minimum zoom
+
+    Parameters
+    ----------
+    precision: string
+        precision level
+
+    Returns
+    -------
+        zoom level
+
+    """
+    if precision == "10m":
+        return 6
+    elif precision == "1m":
+        return 11
+    elif precision == "30cm":
+        return 14
+    else:
+        return 17
+
+
+def _tile2lng(tile_x, zoom):
+    """Returns tile longitude
+
+    Parameters
+    ----------
+    tile_x: int
+        x coordinate
+    zoom: int
+        zoom level
+
+    Returns
+    -------
+        longitude
+    """
+    return ((tile_x / 2 ** zoom) * 360.0) - 180.0
+
+
+def _tile2lat(tile_y, zoom):
+    """Returns tile latitude
+
+    Parameters
+    ----------
+    tile_y: int
+        y coordinate
+    zoom: int
+        zoom level
+
+    Returns
+    -------
+        latitude
+    """
+    n = np.pi - 2 * np.pi * tile_y / 2 ** zoom
+    return (180.0 / np.pi) * np.arctan(0.5 * (np.exp(n) - np.exp(-n)))
+
+
+def _calculate_tile_area(tile):
+    """Returns tile area in square kilometers
+
+    Parameters
+    ----------
+    tile: list
+        tile in format [x,y,z]
+
+    Returns
+    -------
+        area of tile
+
+    """
+    EARTH_RADIUS = 6371.0088
+    left = np.deg2rad(_tile2lng(tile[:, 0], tile[:, 2]))
+    top = np.deg2rad(_tile2lat(tile[:, 1], tile[:, 2]))
+    right = np.deg2rad(_tile2lng(tile[:, 0] + 1, tile[:, 2]))
+    bottom = np.deg2rad(_tile2lat(tile[:, 1] + 1, tile[:, 2]))
+    return (
+        (np.pi / np.deg2rad(180))
+        * EARTH_RADIUS ** 2
+        * np.abs(np.sin(top) - np.sin(bottom))
+        * np.abs(left - right)
+    )
+
+
+def calculate_tiles_area(features, precision):
+    """Calculates the area of tiles
+
+    Parameters
+    ----------
+    features: list
+        features from GeoJSON sources and coordinates
+    precision: string
+        precision level
+
+    Returns
+    -------
+        total area of all tiles in square kilometers
+    """
+    zoom = _convert_precision_to_zoom(precision)
+    tiles = burn(features, zoom)
+    return np.sum(_calculate_tile_area(tiles))

requirements.txt

@@ -5,3 +5,5 @@ requests==2.21.0
 requests-toolbelt==0.9.1
 jsonschema==3.0.1
 jsonseq==1.0.0
+mercantile==1.1.6
+supermercado==0.2.0

setup.py

@@ -34,6 +34,8 @@ def read(fname):
         "requests-toolbelt",
         "jsonschema~=3.0",
         "jsonseq~=1.0",
+        "mercantile~=1.1.6",
+        "supermercado~=0.2.0",
     ],
     include_package_data=True,
     zip_safe=False,