Creating SphericalPolyginsSet for FootprintOverlapDetector

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!

Hi @lghizoni

Welcome to the Orekit

I’m not familiar with that Orekit functionalities. However, I can try to help you.

First of all, I think that discussions in the post below can help you a lot. Especially Manny’s remark about the definition of the points of the polygon is very important (i.e. in a CCW direction).

The Orekit user in this post encountered the same problem that you currently have. If you don’t find any improvement, please find some remarks below.

First, could you change the value of hyperplaneThickness parameter from 10.0 to 1.0e-10. It is the typical value to use.

Secondly, maybe I’m wrong but your are calculated phi angle of S2Point using radians(lat + 90). Could you try using radians(90 - lat) ?

Finally, if my two previous remarks don’t help you, could you try to initialize the SphericalPolygonsSet using the static methods in EllipsoidTessellator class?

geodeticPoints  = []
for lon, lat in zip(longitude, latitude):
  geodeticPoints.append(GeodeticPoint(radians(lat), radians(lon), 0.0)

hyperplaneThickness = 1.0e-10
region_polygon = EllipsoidTessellator.buildSimpleZone(hyperplaneThickness, geodeticPoints)

I hope that I will help you. If not, we will try to find another solution.

Best regards,
Bryan

1 Like

Dear @bcazabonne ,

Thank you so much for the inputs!

I will try implementing your suggestions this evening, and will update you with the results.

Cheers!

Best regards,
/Leonardo.

Dear @bcazabonne ,

I have tried implementing your suggestions to my code. However, I still get weird error messages.

Here’s what I tried, all separately:

  1. Changed hyperplaneThickness to 1.0e-10;
  2. Changed phi from “lat+90” to “90-lat”;
  3. The above two changes together;
  4. Changed my approach to build the region_polygon to the one you suggested, EllipsoidTessellator.buildSimpleZone, together with hyperplaneThickness to 1.0e-10;

All the above changes generate an OutOfMemoryError:

Traceback (most recent call last):
  File "fov_country_event.py", line 194, in <module>
    countryDetector = FootprintOverlapDetector(imagerFOV, earth, region_polygon, sampling_step)
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
java.lang.OutOfMemoryError: GC overhead limit exceeded
        at org.orekit.bodies.Ellipse.projectToEllipse(Ellipse.java:261)
        at org.orekit.bodies.Ellipse.getCenterOfCurvature(Ellipse.java:341)
        at org.orekit.models.earth.tessellation.Mesh$Node.move(Mesh.java:623)
        at org.orekit.models.earth.tessellation.Mesh.addNode(Mesh.java:241)
        at org.orekit.models.earth.tessellation.EllipsoidTessellator.addNode(EllipsoidTessellator.java:594)
        at org.orekit.models.earth.tessellation.EllipsoidTessellator.addAllNeighborsIfNeeded(EllipsoidTessellator.java:582)
        at org.orekit.models.earth.tessellation.EllipsoidTessellator.neighborExpandMesh(EllipsoidTessellator.java:332)
        at org.orekit.models.earth.tessellation.EllipsoidTessellator.sample(EllipsoidTessellator.java:249)
        at org.orekit.propagation.events.FootprintOverlapDetector.sample(FootprintOverlapDetector.java:200)
        at org.orekit.propagation.events.FootprintOverlapDetector.<init>(FootprintOverlapDetector.java:127)

After reading more carefully the documentation, I see that the polygon needs to be defined counter-clock wise, so I simply reversed the latitude and longitude lists as an additional try:

latitude.reverse()
longitude.reverse()

Still, the same kind of error OutOfMemoryError pops out.

Finally, I noticed that a fundamental parameter was too small, the samplingStep for the detector!!

sampling_step = 10.0 # [meters]
countryDetector = FootprintOverlapDetector(imagerFOV, earth, region_polygon, sampling_step)

Changing it to 1000.0, or even 500.0, solved the issue! No more errors :slight_smile:

Thank you so much again for your help! I am learning out how rich and powerful Orekit is, and your suggestion for building the region_polygon with the EllipsoidTessellator.buildSimpleZone just made my life much easier!

I am attaching my code here in case it helps anyone in the future.
fov_country_event.py (6.3 KB)

Best wishes,
/Leonardo.

1 Like

How many points are in your country boundary definition?
SphericalPolygons are not designed to support hundreds of points. In many cases, it
is sufficient to smooth the boundary by decimating the list to something realistic (say less than 100 points).

Also could you tell what are the parameters for tiling? The tile length and width are expected in meters, just like everything in Orekit, so you may try to define too small tiles. Note that the quantization parameter used for fine positioning improves the accuracy but also increases a lot computation. Basically, if you set this parameter to 4 (the recommended value), it means tiles with a length of 60000m and width of 30000m will in fact be computed using points separated by 60000/4=15000m along the tile and 30000/4=7500m accross the track, so the base grid will compute 16 points for each tile, and then select 4 to define the corners, dropping the other 12 points. This means computation and memory cost increases in quantization².

1 Like

Dear @luc ,

Thanks a lot for your inputs.

Currently getting the coordinates for a small country such as Thailand, I have 64 points for the boundary. However, for a bigger country, say Brazil, there are 203 points! Is there a better method to define such big boundaries? You mention smoothing the boundary…is there a method in Orekit to deal with such smoothing? Or should I figure this out by myself, let’s say by getting regional centroids of a sublist of points, and thereby reducing the global amount? Or even by simply drawing by hand a polygon with many less points?

I am not really sure what you mean by parameters for tiling? I’m sorry, as I mentioned before, I am quite a newbie in Orekit.

Thank you once again.

Best regards,
/Leonardo.

Hi @lghizoni

Forget about the tiling quantization, it is more related to tiles than to grids (EllipsoidTessellator manages both aspects, which are dual to each other). The sampling_step is indeed equivalent to tile dimension dived by quantization step, so this was really the culprit as you found by yourself.
Beware that the sampling step is along one dimension and the FootprintOverlapDetector builds a two-dimensional grid to perform its computation, so the number of points in this grid will increase as the inverse of the square of the step size. For a country as large as Brazil, I think 500m would be far too small! You would certainly end up with millions of grid points.

Anyway, FootprintOverlapDetector performs a complex computation and will always be costly. We combine two arbitrary cones centered on two different points (Earth center and spacecraft), its difficult.

I guess 203 points would still be manageable. There are no methods in Orekit or Hipparchus to help smoothing, you’ll have to do it by yourself. The idea is more to avoid tiny details (like thin peninsula or fjords) that would generate noisy boundaries with small enclosed area).

1 Like

Hello @luc

Thanks a lot for the clarifications. I’ve noticed how costly is the FootprintOverlapDetector - I’m having to approximate to some simpler polygons, as I end up getting the same kind of “OutOfMemoryError” when I propagate and try to log the events.

Thanks a lot for your inputs.

/Leonardo.