import spiceypy as spice
import numpy as np
from geojson import Polygon, Feature, FeatureCollection
import geojson
[docs]def lon360to180(lon):
new_lon = lon % 360
new_lon = new_lon - 360 if new_lon > 180 else new_lon
return new_lon
[docs]def roi_pnt2vtx(center_lon, center_lat, diameter, body_radius, n_vertices):
vertex_coords = []
# ROI center point pointing (normalised) vector in body-fixed frame
p_vec = spice.latrec(1.0, lon360to180(center_lon)*spice.rpd(), center_lat*spice.rpd())
# ROI angular diameter
theta = diameter / ( np.cos(center_lat*spice.rpd()) * body_radius )
# fist vertex (rotation of p_vec around +Z axis of body-fixed frame by theta/2 angle)
vertex = spice.vrotv(p_vec, [0,0,1], theta/2)
r, lon, lat = spice.reclat(vertex)
vertex_coords = [ ( lon*spice.dpr(), lat*spice.dpr() ) ]
# define vertices rotation angles around p_vec
alphas = [ (i+1)*(spice.twopi()/n_vertices) for i in range(n_vertices-1)]
# define following vertices
for alpha in alphas:
this_vertex = spice.vrotv(vertex, p_vec, alpha)
r, lon, lat = spice.reclat(this_vertex)
vertex_coords.append( ( lon*spice.dpr(), lat*spice.dpr() ) )
vertex_coords = [ vertex_coords + [vertex_coords[0]] ]
return vertex_coords
[docs]def nomenclature_to_roi(input_geojson_file, output_geojson_file, body_radius, n_vertices):
"""
Turn an input nomenclature GeoJSON file containing the longitude, latitude, and diameter of each POINT feature into an output GeoJSON file containing a 'polygonized' geometry features.
"""
with open(input_geojson_file, 'r') as f:
input_geojson = json.load(f)
ids = []
center_lons = []
center_lats = []
diameters = []
for feature in input_geojson['features']:
ids.append(feature['properties']['name'])
center_lons.append(feature['properties']['center_lon'])
center_lats.append(feature['properties']['center_lat'])
diameters.append(feature['properties']['diameter'])
collection = []
for id, lon, lat, diameter in zip(ids, center_lons, center_lats, diameters):
vertex_coords = roi_pnt2vtx(lon, lat, diameter, body_radius, n_vertices)
feature = Feature(geometry=Polygon(vertex_coords), properties={'id':id, 'diameter':diameter})
collection.append(feature)
with open(output_geojson_file, 'w') as f:
geojson.dump(FeatureCollection(collection), f)
[docs]def add_roi_to_geojson(geojson_file, id, center_lon, center_lat, diameter, body_radius, n_vertices):
with open(geojson_file, 'r') as f:
rois_geojson = geojson.load(f)
# new ROI feature
vertex_coords = roi_pnt2vtx(center_lon, center_lat, diameter, body_radius, n_vertices)
feature = Feature(geometry=Polygon(vertex_coords), properties={'id':id, 'diameter':diameter})
features = rois_geojson['features']
features.append(feature)
new_rois_geojson = FeatureCollection(features)
with open(geojson_file, 'w') as fout:
geojson.dump(new_rois_geojson, fout)
# -- Polygon split implementation by PawaritL
# https://gist.github.com/PawaritL/ec7136c0b718ca65db6df1c33fd1bb11
# https://towardsdatascience.com/around-the-world-in-80-lines-crossing-the-antimeridian-with-python-and-shapely-c87c9b6e1513
from typing import Union, List, Generator
from shapely.geometry import mapping, Polygon, GeometryCollection
from shapely import affinity
[docs]def check_crossing(lon1: float, lon2: float, validate: bool = False, dlon_threshold: float = 180.0):
"""
Assuming a minimum travel distance between two provided longitude coordinates,
checks if the 180th meridian (antimeridian) is crossed.
"""
if validate and any([abs(x) > 180.0 for x in [lon1, lon2]]):
raise ValueError("longitudes must be in degrees [-180.0, 180.0]")
return abs(lon2 - lon1) > dlon_threshold
[docs]def translate_polygons(geometry_collection: GeometryCollection,
output_format: str = "geojson") -> Generator[
Union[List[dict], List[Polygon]], None, None
]:
for polygon in geometry_collection.geoms: # added .geoms
(minx, _, maxx, _) = polygon.bounds
if minx < -180:
geo_polygon = affinity.translate(polygon, xoff=360)
elif maxx > 180:
geo_polygon = affinity.translate(polygon, xoff=-360)
else:
geo_polygon = polygon
yield_geojson = output_format == "geojson"
yield json.dumps(mapping(geo_polygon)) if (yield_geojson) else geo_polygon
import math
import copy
import json
from typing import Union, List
from shapely.geometry import Polygon, LineString, GeometryCollection
from shapely.ops import split
[docs]def split_polygon(geojson: dict, output_format: str = "geojson", validate: bool = False) -> Union[
List[dict], List[Polygon], GeometryCollection
]:
"""
Given a GeoJSON representation of a Polygon, returns a collection of
'antimeridian-safe' constituent polygons split at the 180th meridian,
ensuring compliance with GeoJSON standards (https://tools.ietf.org/html/rfc7946#section-3.1.9)
Assumptions:
- Any two consecutive points with over 180 degrees difference in
longitude are assumed to cross the antimeridian
- The polygon spans less than 360 degrees in longitude (i.e. does not wrap around the globe)
- However, the polygon may cross the antimeridian on multiple occasions
Parameters:
geojson (dict): GeoJSON of input polygon to be split. For example:
{
"type": "Polygon",
"coordinates": [
[
[179.0, 0.0], [-179.0, 0.0], [-179.0, 1.0],
[179.0, 1.0], [179.0, 0.0]
]
]
}
output_format (str): Available options: "geojson", "polygons", "geometrycollection"
If "geometrycollection" returns a Shapely GeometryCollection.
Otherwise, returns a list of either GeoJSONs or Shapely Polygons
validate (bool): Checks if all longitudes are within [-180.0, 180.0]
Returns:
List[dict]/List[Polygon]/GeometryCollection: antimeridian-safe polygon(s)
"""
output_format = output_format.replace("-", "").strip().lower()
coords_shift = copy.deepcopy(geojson["coordinates"])
shell_minx = shell_maxx = None
split_meridians = set()
splitter = None
for ring_index, ring in enumerate(coords_shift):
if len(ring) < 1:
continue
else:
ring_minx = ring_maxx = ring[0][0]
crossings = 0
for coord_index, (lon, _) in enumerate(ring[1:], start=1):
lon_prev = ring[coord_index - 1][0] # [0] corresponds to longitude coordinate
if check_crossing(lon, lon_prev, validate=validate):
direction = math.copysign(1, lon - lon_prev)
coords_shift[ring_index][coord_index][0] = lon - (direction * 360.0)
crossings += 1
x_shift = coords_shift[ring_index][coord_index][0]
if x_shift < ring_minx: ring_minx = x_shift
if x_shift > ring_maxx: ring_maxx = x_shift
# Ensure that any holes remain contained within the (translated) outer shell
if (ring_index == 0): # by GeoJSON definition, first ring is the outer shell
shell_minx, shell_maxx = (ring_minx, ring_maxx)
elif (ring_minx < shell_minx):
ring_shift = [[x + 360, y] for (x, y) in coords_shift[ring_index]]
coords_shift[ring_index] = ring_shift
ring_minx, ring_maxx = (x + 360 for x in (ring_minx, ring_maxx))
elif (ring_maxx > shell_maxx):
ring_shift = [[x - 360, y] for (x, y) in coords_shift[ring_index]]
coords_shift[ring_index] = ring_shift
ring_minx, ring_maxx = (x - 360 for x in (ring_minx, ring_maxx))
if crossings: # keep track of meridians to split on
if ring_minx < -180: split_meridians.add(-180)
if ring_maxx > 180: split_meridians.add(180)
n_splits = len(split_meridians)
if n_splits > 1:
raise NotImplementedError(
"""Splitting a Polygon by multiple meridians (MultiLineString)
not supported by Shapely"""
)
elif n_splits == 1:
split_lon = next(iter(split_meridians))
meridian = [[split_lon, -90.0], [split_lon, 90.0]]
splitter = LineString(meridian)
shell, *holes = coords_shift if splitter else geojson["coordinates"]
if splitter:
split_polygons = split(Polygon(shell, holes), splitter)
else:
split_polygons = GeometryCollection([Polygon(shell, holes)])
geo_polygons = list(translate_polygons(split_polygons, output_format))
if output_format == "geometrycollection":
return GeometryCollection(geo_polygons)
else:
return geo_polygons
import os
from gfinder.geometry import GeometryFactory
[docs]def autoclass_content(module_name, class_name):
content = class_name + '\n'
content += '^' * len(class_name) + '\n'
content += '\n'
content += f'.. autoclass:: gfinder.{module_name}.{class_name}\n'
content += ' :members:\n'
content += ' :undoc-members:\n'
content += ' :show-inheritance:\n'
content += '\n'
return content
[docs]def gen_apidocs(dirpath):
modules = {
'config': 'Config',
'datastore': 'DataStore',
'dataviewer': 'DataViewer',
'detector': 'Detector',
'event': 'Event',
'geometry': 'Geometry',
'geometryevent': 'GeometryEvent',
'opportunity': 'Opportunity',
'radiometry': 'Radiometry',
'resource': 'Resource',
'roi': 'ROI',
'scenario': 'Scenario',
'simulator': 'Simulator',
'tamn': 'TAMN',
'target': 'Target',
'timeline': 'Timeline',
'utils': 'Utils' }
geometry_factory = GeometryFactory()
core_geometry_classes = ['Geometry', 'Condition', 'GeometryFactory']
lowlevel_class_names = geometry_factory.lowlevel_class_names
highlevel_class_names = geometry_factory.highlevel_class_names
excluded_class_names = core_geometry_classes + lowlevel_class_names + highlevel_class_names
all_geometries = geometry_factory.get_geometries()
all_geometries.sort()
for key in modules.keys():
title = f'``{modules[key]}`` module'
content = f'.. _ref_gfinder_{key}_module:\n'
content += '\n'
content += f'{title}\n'
content += '=' * len(title) + '\n'
content += '\n'
if key == 'geometry':
content += 'Core classes\n'
content += '------------\n'
content += '\n'
for class_name in core_geometry_classes:
content += autoclass_content(key, class_name)
content += 'Low-level Geometry classes\n'
content += '--------------------------\n'
content += '\n'
for class_name in geometry_factory.lowlevel_class_names:
content += autoclass_content(key, class_name)
content += 'High-level Geometry classes\n'
content += '---------------------------\n'
content += '\n'
for class_name in geometry_factory.highlevel_class_names:
content += autoclass_content(key, class_name)
content += 'Geometry classes\n'
content += '----------------\n'
content += '\n'
for class_name in all_geometries:
if class_name not in excluded_class_names:
content += autoclass_content(key, class_name)
else:
content += f'.. automodule:: gfinder.{key}\n'
content += ' :members:\n'
content += ' :undoc-members:\n'
content += ' :show-inheritance:'
# delete rst file if exists
rstfile = f'{dirpath}/{key}.rst'
if os.path.exists(rstfile):
os.remove(rstfile)
print(f'removed {rstfile} file.')
# write rst file
with open(rstfile, 'w') as f:
f.write(content)
print(f'written {rstfile} file.')
