Extrapolate GeographicZoneDetector results

Orekit Python Wrapper usage
1

Hello, i am new here and new in Python world.
I am writing a code to Print the coordinates/time of the spacecraft when enter/leave a given country. Everythig works: i have my satellites with the GroundTrack coordinates in a period of 12 hours (thanks get_pv_coordinates from TLEs), i have the polygons representing the country (thanks EllipsoidTessellator) , i also initialized GeographicZoneDetector (countryDetector = GeographicZoneDetector(earth, region_polygon, margin)).
The problem is that i don’t know how to show the result of the GeographicZoneDetector (coordinates/time of the spacecraft when enter/leave a given country)
Can anyone help me in this little task ? Please

Thanks a lot.

2

Hi @rekoraffa welcome,

You should set up an Earth shape model (with OneAxisEllipsoid) which I guess you already have as it is needed for GeographicZoneDetector.
Then you should get the SpacecraftState at event occurrence. You either get it if you implement the eventOccurred method in a custom EventHandler or you get it after the fact if you logged all events using EventLogger.

Then, from the Earth model and the spacecraft state, you can build a GeodeticPoint using (in Java):

  GeodeticPoint gp = earth.transform(state.getPVCoordinates().getPosition(),
                                     state.getFrame(),
                                     state.getDate());
3

Thanks a lot to spend time for me, i really appreciate.
Regarding the SpacecraftState, i have the PV vector propaganted from TLEs. It is enought to determin SpacecraftState

4

You can get the full spacecraft state already set up by the TLE propagator either in the eventOccurred method or if you dumped that and only kept the date you can just call again tlepropagator.propagate(date) to retrieve it again.

5

This is my code. it works (no errors occur) but i don’t know why i am not able to show the result of GeographicZoneDetector. can you help me ?
Thanks

utc = TimeScalesFactory.getUTC()
t0  = AbsoluteDate(2020,  8,  1, 21, 0, 0.0, utc)
t1  = AbsoluteDate(2020,  8,  2, 4, 0, 0.0, utc)

#import TLE 
tles = TLESeries.from_file(filename)

# Set the reference frame where the satellite is defined.
j2000 = FramesFactory.getEME2000()
inertialFrame = j2000 

#For a given interval time and TLE series, POSITIONS and VELOCITIES are extrapolated within a the specific reference frame.
pv = tles.get_pv_coordinates(t0=t0, t1=t1, step = stepTime , frame=j2000)

#CALCULATE GEOGRAPHICAL SUBPOINT

# Overall duration for extrapolation
startDate = t0
duration = 12 * 60 * 60
# Set the step time 
step_time = 20.
# Time array in orekit AbsoluteDate format
t = [startDate.shiftedBy(float(dt)) \
        for dt in np.arange(0, duration, step_time)]

#Create a planet to map positions on
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                         Constants.WGS84_EARTH_FLATTENING, 
                         ITRF)

#Create the satellite position and velocity vector (3-D vector) from pv
p = [tpv.getPosition() for tpv in pv]

#Calculate the geographical subpoint. Each subpoints requires a position from the p-vector and a time stamp from the t-vector
subpoint = [earth.transform(tp, j2000, tt)  for tt, tp in zip(t,p)]  #Cartesian point
#From geographical subpoint extrapolate LON, LAT, ALT
lat = np.degrees([gp.getLatitude()  for gp in subpoint])
lon = np.degrees([gp.getLongitude() for gp in subpoint])

# Region-Of-Interest Specification
region_of_interest = 'Libya'

#USING CARTYOPY: Get Region-Of-Interest Polygon - Available Resolutions: 110m, 50m, and 10m
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

# List of Geographical Coordinates of the Polygon vertices
longitude, latitude = country_geometry.exterior.coords.xy

# Prepare latitude and longitude to specify Region-Of-Interest Polygon
latitude.reverse()
longitude.reverse()

# Create Geodetic Points
geodeticPoints = []
for lon, lat in zip(longitude, latitude):
	#print('{}   {}'.format(lat, lon))
	geodeticPoints.append(GeodeticPoint(radians(lat), radians(lon), 0.0))

# Build Region Polygon
hyperplaneThickness = 1.0e-10
region_polygon = EllipsoidTessellator.buildSimpleZone(hyperplaneThickness, geodeticPoints)

#USING OREKIT DETECTOR
margin= 0.0     #Angular margin to apply to the zone
countryDetector = GeographicZoneDetector(earth, region_polygon, margin)  #GeographicZoneDetector

ZoneDetectorLogger = EventsLogger()
logged_detector = ZoneDetectorLogger.monitorDetector(countryDetector)

# We add the Zone Detector, together with the eventslogger to our propagator.
propagator = KeplerianOrbit(pv, j2000, t0, Constants.WGS84_EARTH_MU)
propagator.addEventDetector(logged_detector)

start_time = None
contacts = pd.DataFrame()
events = ZoneDetectorLogger.getLoggedEvents()
for event in ZoneDetectorLogger.getLoggedEvents():
    if event.isIncreasing():
        start_time = event.getState().getDate()
    elif start_time:
        stop_time = event.getState().getDate()
        contacts = contacts.append({"Duration (s)": stop_time.durationFrom(start_time),
                                    "Start Time (UTC)": absolutedate_to_datetime(start_time),
                                    "Stop Time (UTC)": absolutedate_to_datetime(stop_time)},
                                    ignore_index=True)
        contacts.to_csv('Accesses', sep='\t', index=False)
        start_time = None
print(contacts)