Skip to content

Commit

Permalink
Add GeoJSON Area Calc (#19)
Browse files Browse the repository at this point in the history
  • Loading branch information
nllong authored Feb 21, 2023
1 parent ae879c9 commit df67582
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 6 deletions.
7 changes: 4 additions & 3 deletions pyseed/seed_client.py
Original file line number Diff line number Diff line change
Expand Up @@ -229,10 +229,11 @@ def get_buildings(self) -> List[dict]:
return buildings

def get_property_view(self, property_view_id: int) -> dict:
"""Return a single property (view and state) by the property view id.
"""Return a single property (view and state) by the property view id. It is
recommended to use the more verbose version of `get_property` below.
Args:
property_view_id (int): ID of the property to return. This is the ID that is in the URL http://SEED_URL/app/#/properties/{property_view_id}
property_view_id (int): ID of the property to return. This is the ID that is in the URL http://SEED_URL/app/#/properties/{property_view_id} and resolves to {host}/api/v3/property_views/{property_view_id}
Returns:
dict: {
Expand Down Expand Up @@ -440,7 +441,7 @@ def update_labels_of_buildings(
Args:
add_label_names (list): list of label names to add, will be converted to IDs
remove_label_names (list): list of label names to remove, will be converted to IDs
building_ids (list): list of building IDs to add/remove labels
building_ids (list): list of building IDs (property_view_id) to add/remove labels
inventory_type (str, optional): taxlot or property inventory. Defaults to 'property'.
Raises:
Expand Down
97 changes: 97 additions & 0 deletions pyseed/utils.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,104 @@
# from __future__ import division

# Imports from Third Party Modules
import csv
import json
from math import pi, sin
from pathlib import Path

WGS84_RADIUS = 6378137


def _rad(value):
return value * pi / 180


def _ring_area(coordinates):
"""
Calculate the approximate total_area of the polygon were it projected onto
the earth. Note that this _area will be positive if ring is oriented
clockwise, otherwise it will be negative.
Reference:
Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for
Polygons on a Sphere", JPL Publication 07-03, Jet Propulsion
Laboratory, Pasadena, CA, June 2007 http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
@Returns
{float} The approximate signed geodesic total_area of the polygon in square meters.
"""

assert isinstance(coordinates, (list, tuple))

total_area = 0
coordinates_length = len(coordinates)

if coordinates_length > 2:
for i in range(0, coordinates_length):
if i == (coordinates_length - 2):
lower_index = coordinates_length - 2
middle_index = coordinates_length - 1
upper_index = 0
elif i == (coordinates_length - 1):
lower_index = coordinates_length - 1
middle_index = 0
upper_index = 1
else:
lower_index = i
middle_index = i + 1
upper_index = i + 2

p1 = coordinates[lower_index]
p2 = coordinates[middle_index]
p3 = coordinates[upper_index]

total_area += (_rad(p3[0]) - _rad(p1[0])) * sin(_rad(p2[1]))

total_area = total_area * WGS84_RADIUS * WGS84_RADIUS / 2

return total_area


def _polygon_area(coordinates):

assert isinstance(coordinates, (list, tuple))

total_area = 0
if len(coordinates) > 0:
total_area += abs(_ring_area(coordinates[0]))

for i in range(1, len(coordinates)):
total_area -= abs(_ring_area(coordinates[i]))

return total_area


def geojson_area(geometry):
"""Calculate the area of a GeoJSON feature. This method is taken from
a combination of ChatGPT conversion of:
https://github.com/mapbox/geojson-area/blob/master/index.js
and
https://github.com/scisco/area/blob/master/area/__init__.py"""

if isinstance(geometry, str):
geometry = json.loads(geometry)

assert isinstance(geometry, dict)

total_area = 0

if geometry['type'] == 'Polygon':
return _polygon_area(geometry['coordinates'])
elif geometry['type'] == 'MultiPolygon':
for i in range(0, len(geometry['coordinates'])):
total_area += _polygon_area(geometry['coordinates'][i])
elif geometry['type'] == 'GeometryCollection':
for i in range(0, len(geometry['geometries'])):
total_area += geojson_area(geometry['geometries'][i])

return total_area


def read_map_file(mapfile_path):
"""Read in the mapping file"""
Expand Down
7 changes: 4 additions & 3 deletions tests/test_seed_client.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,11 +103,12 @@ def test_search_buildings(self):
assert prop["state"]["address_line_1"] == "111 Street Lane, Chicago, IL"
assert prop["state"]["extra_data"]["EUI Target"] == 120

# test the property view (same as previous, just less data)
# test the property view (same as previous, just less data). It
# is recommended to use `get_property` instead.
prop = self.seed_client.get_property_view(properties[0]["id"])
print(prop)
assert prop["id"] == properties[0]["id"]
assert prop["state"]["address_line_1"] == "111 Street Lane, Chicago, IL"
assert prop["state"]["extra_data"]["EUI Target"] == 120
assert prop["cycle"]["name"] == "pyseed-api-test"

# There are 2 if filtering, because B-1 and B-10
properties = self.seed_client.search_buildings(identifier_filter="B-1")
Expand Down

0 comments on commit df67582

Please sign in to comment.