Hello All,
I am trying to create a SphericalPolygonsSet to be inserted in a FootprintOverlapDetector in python. Basically, I get latitude and longitude coordinates of a given country from cartopy, transform them to theta and phi,create a list of S2Points, and insert it to the SphericalPolygonsSet. So far so good…the thing I’m not being able to handle with is the hyperplaneThickness, parameter which I’m quite confused on how to specify. I see from the documentation that a rule of thumb is " a value slightly below the size of the most accurate detail needed is a good value for the hyperplaneThickness parameter". So tried values between 100.0 and 1.0, and even very small ones such as 1.0e-3, but always getting the following error (plus, I’m getting errors for geographical coordinates that shouldn’t be there, since I’m dealing with the Thailand territory):
File “fov_country_event.py”, line 171, in
countryDetector = FootprintOverlapDetector(imagerFOV, earth, region_polygon, sampling_step)
orekit.JavaError: <super: <class ‘JavaError’>, >
Java stacktrace:
org.orekit.errors.OrekitException: cannot compute aiming direction at singular point: latitude = -47.762, longitude = -31.49
at org.orekit.models.earth.tessellation.Mesh.(Mesh.java:104)
at org.orekit.models.earth.tessellation.EllipsoidTessellator.sample(EllipsoidTessellator.java:243)
at org.orekit.propagation.events.FootprintOverlapDetector.sample(FootprintOverlapDetector.java:200)
at org.orekit.propagation.events.FootprintOverlapDetector.(FootprintOverlapDetector.java:127)
Here’s my code so far:
import orekit
vm = orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
from org.orekit.time import TimeScalesFactory
from org.orekit.time import AbsoluteDate
from org.orekit.utils import Constants
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.propagation.analytical import EcksteinHechlerPropagator
from org.orekit.orbits import OrbitType
from org.orekit.propagation.analytical import KeplerianPropagator
from org.hipparchus.geometry.euclidean.threed import Line
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.frames import FramesFactory
from org.orekit.utils import IERSConventions
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.geometry.fov import DoubleDihedraFieldOfView
from org.hipparchus.geometry.spherical.twod import S2Point
from org.hipparchus.geometry.spherical.twod import SphericalPolygonsSet
from org.orekit.propagation.events import EventDetector, EventsLogger
from org.orekit.propagation.events import FootprintOverlapDetector
from org.orekit.attitudes import NadirPointing
import math
from math import radians, degrees, atan, pi
from IPython.display import FileLink, display
import matplotlib.pyplot as plt
import datetime
from beyond.dates import Date
from beyond.utils.ltan import ltan2raan
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.shapereader as shpreader
import numpy as np
import pandas as pd
import sys
sys.path.append('../ORB/')
import orbit_properties as op
R_earth = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
Mu_earth = Constants.WGS84_EARTH_MU
utc = TimeScalesFactory.getUTC()
begin_date = datetime.datetime(2024, 5, 1, 10, 30, 0) # (YYYY,M,D,H,M,S)
date = Date(begin_date.year,begin_date.month,begin_date.day,begin_date.hour,begin_date.minute,begin_date.second)
epochDate = AbsoluteDate(begin_date.year,begin_date.month,begin_date.day,begin_date.hour,begin_date.minute,float(begin_date.second),utc)
initialDate = epochDate
duration = 2 # Days
step = 10.0
ra = 500 * 1000 # Apogee
rp = 500 * 1000 # Perigee
a = (ra + rp + 2*Constants.WGS84_EARTH_EQUATORIAL_RADIUS)/2.0 # Semi-Major Axis
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS)/a
i = radians(op.sso_inclination(ra*1e-3,e))
raan = float(ltan2raan(date,(begin_date.hour+begin_date.minute/60),type='mean'))
omega = radians(0.0)
lv = radians(0.0)
swath = 15.0 # km
xfov = degrees(atan((swath/2)/(ra*1e-3)))
yfov = 30.0 + xfov
fov_margin = 0.0
region_of_interest = 'Thailand'
region_colour = 'r'
inertialFrame = FramesFactory.getEME2000()
initialOrbit = KeplerianOrbit(a,e,i,omega,raan,lv,PositionAngle.TRUE,inertialFrame,initialDate,Constants.WGS84_EARTH_MU)
propagator = EcksteinHechlerPropagator(initialOrbit,
Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20,
Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40,
Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60)
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ITRF)
attitudeProvider = NadirPointing(inertialFrame, earth)
propagator.setAttitudeProvider(attitudeProvider)
initialState = propagator.getInitialState()
resolution = 110
shpfilename = shpreader.natural_earth(resolution=str(resolution)+'m',category='cultural', name='admin_0_countries')
reader = shpreader.Reader(shpfilename)
countries = reader.records()
country_geometry = pd.DataFrame()
for country in countries:
if country.attributes['NAME'] == region_of_interest:
country_geometry = country.geometry
longitude, latitude = country_geometry.exterior.coords.xy
theta = []
phi = []
for lon, lat in zip(longitude, latitude):
theta.append(radians(lon))
phi.append(radians(lat + 90))
polygon_vertices = []
for theta_p, phi_p in zip(theta, phi):
polygon_vertices.append(S2Point(theta_p, phi_p))
hyperplaneThickness = 10.0
region_polygon = SphericalPolygonsSet(hyperplaneThickness, polygon_vertices)
imagerFOV = DoubleDihedraFieldOfView(Vector3D.PLUS_K, Vector3D.PLUS_I, radians(xfov), Vector3D.PLUS_J, radians(yfov),fov_margin)
sampling_step = 10.
countryDetector = FootprintOverlapDetector(imagerFOV, earth, region_polygon, sampling_step)
My apologies if this is something trivial. It’s the first time I get hands on OREKIT.
Thank you very much in advance!
