"""TAMN module.
This module is a wrapper to the TAMN web service, and provides additional handling for LineString geometries.
"""
import json
from geojson import LineString, MultiLineString
import spiceypy as spice
import numpy as np
import requests
import geodaisy.converters as convert
TAMN_URL = 'http://tamn.spacefrog.design/'
DEFAULT_RADIUS = 50000 # 50000
"""Default radius used in TAMN algorithm.
"""
DEFAULT_POLE_DISTANCE = 0 # 500000
"""Default pole distance used in TAMN algorithm.
"""
[docs]def split_line(line):
coords = line
lines = []
i0 = 0
for i in range(len(coords)-1):
flag, p0 = crossing(coords[i],coords[i+1])
if flag:
if i0 > 0:
line = []
line.append(p0_next)
for coord in coords[i0:i]:
line.append(coord)
else:
line = coords[i0:i]
line.append(p0)
if len(line) > 1: # exclude one-point line (workaround)
lines.append(line)
p0_next = [-p0[0], p0[1]] # first point of the next line
i0 = i + 1 # start index of next line
# add last line
if i0 > 0:
line = []
line.append(p0_next)
for coord in coords[i0:]:
line.append(coord)
lines.append(line)
else:
lines = [line]
return lines
[docs]def crossing(p1, p2): # [lon, lat]
p0 = [0., 0.]
p1_x = np.cos(p1[0]*spice.rpd())
p2_x = np.cos(p2[0]*spice.rpd())
if p1_x < 0 and p2_x < 0:
p1_y = np.sin(p1[0]*spice.rpd())
p2_y = np.sin(p2[0]*spice.rpd())
if (p1_y > 0 and p2_y < 0) or (p1_y < 0 and p2_y > 0):
lat0 = 0.5*(p1[1]+p2[1])
lon0 = p1_y/abs(p1_y) * 180.
p0 = [lon0, lat0]
return True, p0
return False, p0
[docs]def getSplitLineString(feature):
if feature['geometry']['type'] != 'LineString':
return feature
coords = feature['geometry']['coordinates']
lines = split_line(coords)
feature['geometry'] = MultiLineString(lines)
# print(type(feature))
# if not feature.is_valid:
# print('[WARNING] Invalid GeoJSON Feature.')
return feature
[docs]def getSplitFeature(feature):
# Handle LineString geometries independently of TAMN.
#
if feature['geometry']['type'] == 'LineString':
return getSplitLineString(feature)
# Otherwise, let TAMN handle other geometry types.
#
# set TAMN web server end-point
request = TAMN_URL + 'split'
# convert input feature GeoJSON time_inputs_dict in WKT string
wkt = convert.geojson_to_wkt(json.dumps(feature['geometry']))
data = json.dumps({"geometry":wkt, "radius":DEFAULT_RADIUS, "pole_distance":DEFAULT_POLE_DISTANCE})
# send split request to TAMN web service
r = requests.post(request,
data=data,
headers={
'Content-Type':'application/json'
}
)
if r.status_code != 200:
print('[ERROR] POST failed - {code} {text}'.format(code=r.status_code, text=r.text))
return 0
# set new geometry
feature['geometry'] = json.loads(r.text)['geometry']
return feature
[docs]def splitGeoJSON(in_geojson_file, out_geojson_file):
# Load input GeoJSON file that can contain either only one Feature or a FeatureCollection.
with open(in_geojson_file) as f:
geojson_dict = json.load(f)
# If only one Feature, split and update feature geometry in geojson_dict
if geojson_dict['type'] == 'Feature':
splitFeature = getSplitFeature(geojson_dict)
geojson_dict['geometry'] = splitFeature
# If FeatureCollection, split and update each feature
elif geojson_dict['type'] == 'FeatureCollection':
splitFeatures = []
for feature in geojson_dict['features']:
splitFeatures.append(getSplitFeature(feature))
geojson_dict['features'] = splitFeatures
# Write geojson_dict input ouput GeoJSON file
with open(out_geojson_file, 'w') as f:
json.dump(geojson_dict, f)
