pyorbital
pyorbital copied to clipboard
pytroll
Add script to get SNOs for two platforms
...typically Calipso and an Imager satellite like Metop/NOAA
Signed-off-by: Adam Dybbroe [email protected]
- [ ] Closes #xxxx
- [ ] Tests added
- [x] Tests passed
- [x] Passes
flake8 pyorbital - [x] Fully documented
Coverage decreased (-7.3%) to 76.609% when pulling 91f9f7401391e610e157bdaca4412e3aa30a4fba on adybbroe:add-sno-finder-script into 7efe253e20253f3a28501d4f2a4a86efa33805e4 on pytroll:master.
I have added the SNO script, it works on my laptop, I have a bit of unittests (though they are kind of system-tests) and some documentation. I was considering if we should/could merge it and then I could continue doing refactoring and making the tests more like real unittests, etc? What do you think @pnuu @mraspaud @djhoese ?
@adybbroe I tested this with one day's TLE file which contains CALIPSO (29108) and Sentinel-5P (42969).
TLE file
tle-20200102.txt
1 29108U 06016B 20002.14646990 .00000119 00000-0 32995-4 0 9999
2 29108 98.2112 309.9642 0001231 75.1794 284.9543 14.62368728727882
1 29108U 06016B 20002.83069859 +.00000101 +00000-0 +29299-4 0 9992
2 29108 098.2112 310.6444 0001224 075.1808 284.9529 14.62368891727963
1 42969U 17064A 20002.83823431 -.00000020 +00000-0 +11213-4 0 9999
2 42969 098.7199 304.3661 0001110 087.5775 272.5529 14.19547577115111
snos.yaml
station:
longitude: 80
latitude: 70
altitude: 0.03
tle-dirs:
- ../data/tle/
tle-file-format: 'tle-%Y%m%d.txt'
reference_platform:
name: CALIPSO
Terminal run
$python ./get_snos.py -s 20200102 -e 20200103 -t 2 -p Sentinel-5P -c ./snos.yaml
2020-01-02 13:47, 12.5, 2020-01-02 13:46, 38.7, 11508, 80.03, -99.17, 0.6, False
As the intersect time is around 2020-01-02 13:46, I downloaded two data files:
CALIPSO: CAL_LID_L2_VFM-Standard-V4-20.2020-01-02T13-42-15ZN.hdf
Sentinel-5P: S5P_OFFL_L2__NO2____20200102T122416_20200102T140546_11508_01_010302_20200104T053822.nc
However, I plotted the Sentinel-5P swath and CALIPSO lines:
The intersect isn't at 99.17W and 80.03N.
Datetime in both datasets:
CALIPSO:
[datetime.datetime(2020, 1, 2, 13, 42, 14, 837201, tzinfo=<UTC>)
datetime.datetime(2020, 1, 2, 13, 42, 15, 581199, tzinfo=<UTC>)
datetime.datetime(2020, 1, 2, 13, 42, 16, 325200, tzinfo=<UTC>) ...
datetime.datetime(2020, 1, 2, 14, 28, 14, 305199, tzinfo=<UTC>)
datetime.datetime(2020, 1, 2, 14, 28, 15, 49200, tzinfo=<UTC>)
datetime.datetime(2020, 1, 2, 14, 28, 15, 793199, tzinfo=<UTC>)]
Sentinel-5P:
array(['2020-01-02T12:45:50.235000Z', '2020-01-02T12:45:51.075000Z',
'2020-01-02T12:45:51.915000Z', ..., '2020-01-02T13:44:12.957000Z',
'2020-01-02T13:44:13.797000Z', '2020-01-02T13:44:14.637000Z'],
dtype=object)
13:46 is in the time period of the CALIPSO data, but not for the Sentinel-5P data. Note that the next Sentinel-5P file begins after 14:05.
Debug
I tried to output some variables for debugging:
The tobj, delta_t before get_sno_point is 2020-01-02 12:08:00 and 0:04:00.
There's a large difference between 12:08 and the 13:46 of sno.
Here's the (coord_start, coord_end) in the get_arc function:
(-18.652184966983476, 78.69164595153211) (-103.4240079209052, 79.44240885811259)
(-93.21301074810981, 80.50824312832857) (-135.17883694428102, 70.07349226193958)
It seems the first one (ref_pltfrm: CALIPSO) runs so fast. Could that be the problem?
We hade som problems that this script is not accurate enough as the approximation of the path is not exact enought and done with only one arc. I made a more accurate version with several smaller arcs (see below). It takes more time but better finds the SNOs.
# -*- coding: utf-8 -*-
import sys
import os
from datetime import datetime, timedelta
import pyresample as pr
from pyresample.spherical import SCoordinate, Arc
#from pyresample.spherical_geometry import Coordinate as SGCoordinate
#from pyresample.spherical_geometry import Arc as SGArc
#from pyorbital.etc.platforms
from pyorbital.orbital import Orbital
import numpy as np
from pyorbital.tlefile import Tle
import time
t= time.time()
import logging
LOG = logging.getLogger('snos')
handler = logging.StreamHandler(sys.stderr)
handler.setLevel(0)
LOG.setLevel(0)
LOG.addHandler(handler)
# Norrköping coordinates:
NRK_LON = 16.1465
NRK_LAT = 58.5780
NRK_ALT = 0.03
R = 6370997.0
PATTERN = [0,1,-1,2,-2]
#TLE_DIR = '/data/orbital_elements/TLE/'
TLE_DIR = '/home//temp'
TLE_SUBD_FORMAT = '%Y%m'
TLE_FILE_FORMAT = 'tle-%Y%m%d????.txt'
TLE_FILE_FORMAT2 = 'tle-%Y%m%d.txt'
TLE_FILE_FORMAT3 = 'TLE_noaa18.txt'
class NoTleFile(Exception):
pass
TLE_BUFFER_OTHER = {}
TLE_BUFFER_CALIPSO = {}
TLE_SATNAME = {'npp' : 'SUOMI NPP',
'aqua' : 'AQUA',
'metopb': 'METOP-B',
'metopa': 'METOP-A',
'noaa19' : 'NOAA 19',
'noaa18' : 'NOAA 18',
'sentinel3a' : 'SENTINEL-3A',
'sentinel3b' : 'SENTINEL-3B',
'fengyun3d' : 'FENGYUN 3D',
'noaa15' : 'NOAA 15',
'noaa16' : 'NOAA 16'}
TLE_IDs = {'npp' : 'xxxxx',
'aqua' : 'xxxxx',
'calipso': '29108',
'metopb': '38771',
'metopa': '29499',
'noaa19' : '33591',
'noaa18' : '28654',
'sentinel3a' : 'xxxxx',
'sentinel3b' : 'xxxxx',
'fengyun3d' : 'xxxxx',
'noaa15' : 'xxxxx',
'noaa16' : 'xxxxx'}
max_tle_days_diff = 3
def _get_tle_file(timestamp):
# Find a not too old TLE file
import glob
for delta in PATTERN:
dt_obj = timestamp - timedelta(days=delta)
path = os.path.join(TLE_DIR, dt_obj.strftime(TLE_SUBD_FORMAT))
if os.path.isdir(path):
search = os.path.join(path, dt_obj.strftime(TLE_FILE_FORMAT))
flst = glob.glob(search)
if len(flst) != 0:
fname = flst[0]
LOG.info("Found TLE file: '%s'" % fname)
return fname
elif 1==2:
search = os.path.join(path, dt_obj.strftime(TLE_FILE_FORMAT2))
flst = glob.glob(search)
if len(flst) != 0:
fname = flst[0]
LOG.info("Found TLE file: '%s'" % fname)
return fname
raise NoTleFile("Found no TLE file close in time to " +
str(timestamp.strftime(TLE_FILE_FORMAT)))
def get_tle(platform, timestamp=None):
"""Get the tle from file, if not loaded already"""
stamp = platform + timestamp.strftime('-%Y%m%d')
try:
tle = TLE_BUFFER_OTHER[stamp]
except KeyError:
tle = Tle(TLE_SATNAME[platform], _get_tle_file(timestamp))
TLE_BUFFER_OTHER[stamp] = tle
return tle
def get_tle_archive_calipso(timestamp):
return get_tle_archive(timestamp, filename_calipso, TLE_ID_CALIPSO, TLE_BUFFER_CALIPSO)
def get_tle_archive_other(timestamp):
return get_tle_archive(timestamp, filename_other, TLE_ID_OTHER, TLE_BUFFER_OTHER)
def populate_tle_buffer(filename, TLE_ID, MY_TLE_BUFFER):
with open(filename, 'r') as fh_:
tle_data_as_list = fh_.readlines()
for ind in range(0, len(tle_data_as_list) , 2):
if TLE_ID in tle_data_as_list[ind]:
tle = Tle(TLE_ID, line1=tle_data_as_list[ind], line2=tle_data_as_list[ind+1])
#dto = datetime.strptime(tle.epoch, '%Y-%m-%dT%H:%M:%S:%f')
ts = (tle.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
tobj = datetime.utcfromtimestamp(ts)
MY_TLE_BUFFER[tobj] = tle
def get_tle_archive(timestamp, filename, TLE_ID, MY_TLE_BUFFER):
# read tle data if not already in buffer
if len(MY_TLE_BUFFER) == 0:
populate_tle_buffer( filename, TLE_ID, MY_TLE_BUFFER)
for tobj in MY_TLE_BUFFER:
if tobj > timestamp:
deltat = tobj - timestamp
else:
deltat = timestamp - tobj
if np.abs((deltat).days) < 1:
return MY_TLE_BUFFER[tobj]
for delta_days in range(1, max_tle_days_diff + 1, 1):
for tobj in MY_TLE_BUFFER:
if tobj > timestamp:
deltat = tobj - timestamp
else:
deltat = timestamp - tobj
if np.abs((deltat).days) <= delta_days:
print("Did not find TLE for {:s}, Using TLE from {:s}".format(tobj.strftime("%Y%m%d"),
timestamp.strftime("%Y%m%d")))
return MY_TLE_BUFFER[tobj]
print("Did not find TLE for {:s} +/- 3 days")
def get_arc(timeobj, delta_t, sat):
"""Get the arc defining the sub-satellite track from start time 'timeobj'
to start time 'timeobj' plus time step 'delta_t'"""
# Get the start and end positions:
pos_start = sat.get_lonlatalt(timeobj - delta_t)
pos_end = sat.get_lonlatalt(timeobj + delta_t)
coord_start = pr.spherical.SCoordinate(lon=np.deg2rad(pos_start[0]),
lat=np.deg2rad(pos_start[1]))
coord_end = pr.spherical.SCoordinate(lon=np.deg2rad(pos_end[0]),
lat = np.deg2rad(pos_end[1]))
return pr.spherical.Arc(coord_start, coord_end)
def get_arc_vector(timeobj, delta_t, sat, arc_len_min):
"""Get the arc defining the sub-satellite track from start time 'timeobj'
to start time 'timeobj' plus time step 'delta_t'"""
# Get positions at several points between +/- delta_t:
len_one_arc_s = 60 * arc_len_min # s
num = int(delta_t.seconds/len_one_arc_s)
#print(2*num)
pos_vec = [sat.get_lonlatalt(timeobj + ind * 1.0/num * delta_t) for ind in range(-num, num +1,1)]
# Calculate the Arc for each pixel. Later we sill see if arcs cross each other.
# We could use only start and end. But then the SNO would be approximated some times to much.
arcs = [ pr.spherical.Arc(pr.spherical.SCoordinate(lon=np.deg2rad(point1[0]),
lat=np.deg2rad(point1[1])),
pr.spherical.SCoordinate(lon=np.deg2rad(point2[0]),
lat = np.deg2rad(point2[1]))) for point1, point2 in zip(pos_vec[0:-1], pos_vec[1:])]
return arcs
if __name__ == "__main__":
if len(sys.argv) < 6:
print("Usage: %s " % sys.argv[0] +
("<start-time (yyyymmdd)> <end-time (yyyymmdd)>" +
" <time-window (minutes)> <platform> <calipso platform> <approximation_arc_minutes 2-good 5-fast>"))
sys.exit()
else:
time_start = datetime.strptime(sys.argv[1], "%Y%m%d%H%M")
time_end = datetime.strptime(sys.argv[2], "%Y%m%d%H%M")
minthr = float(sys.argv[3])
platform_id = sys.argv[4]
calipso_id = sys.argv[5]
arc_len_min = int(sys.argv[6])
if platform_id not in TLE_SATNAME:
LOG.error("Platform %s not supported!" % platform_id)
sys.exit()
filename_calipso = "/home/sm_kgkar/Projects/Get_SNOs/TLEs/TLE_{:s}.txt".format(calipso_id)
TLE_ID_CALIPSO = TLE_IDs[calipso_id]
filename_other = "/home/sm_kgkar/Projects/Get_SNOs/TLEs/TLE_{:s}.txt".format(platform_id)
TLE_ID_OTHER = TLE_IDs[platform_id]
import math
minthr_step = 20 # min less than half an orbit probably
dtime = timedelta(seconds = 60 * minthr_step * 2.0)
#timestep_half = timedelta(seconds = 60 * minthr_step * 0.5)
timestep_double = timedelta(seconds = 60 * minthr_step * 2.0)
#timestep = timedelta(seconds = 60 * minthr_step * 1.0)
timestep_plus_30s = timedelta(seconds = 60 * minthr_step * 1.0 + (minthr*0.5)*60 +30) # make sure the two sat pass the SNO in the same step. We need and overlap of at least half minthr minutes.
#delta_t = timedelta(seconds = 1200)
#delta_t = timedelta(seconds = 60 * (minthr_step-1) * 2.0)
tobj = time_start
tle_calipso = None
tle_the_other_one = None
tobj_tmp = time_start
isNorrk = False
i=0
t_diff = timedelta(days=1)
while tobj < time_end:
#print(tobj)
i=i+1
if i==100:
message = time.time()-t, "seconds", dtime
print (message)
if not tle_calipso or calipso_obj - tobj > t_diff or tobj - calipso_obj > t_diff:
tle_calipso = get_tle_archive_calipso(tobj)
if tle_calipso:
ts = (tle_calipso.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
calipso_obj = datetime.utcfromtimestamp(ts)
if not tle_the_other_one or other_obj - tobj > t_diff or tobj - other_obj > t_diff:
tle_the_other_one = get_tle_archive_other(tobj)
if tle_the_other_one:
ts = (tle_the_other_one.epoch - np.datetime64('1970-01-01T00:00:00Z')) / np.timedelta64(1, 's')
other_obj = datetime.utcfromtimestamp(ts)
calipso = Orbital(calipso_id,
line1=tle_calipso.line1,
line2=tle_calipso.line2)
the_other_one = Orbital(platform_id,
line1=tle_the_other_one.line1,
line2=tle_the_other_one.line2)
# Approximated intersection. Track approximated with one arc!
# This is not acurate enough to find all SNOs.
# arc_calipso = get_arc(tobj, timestep_plus_30s, calipso)
# arc_the_other_one = get_arc(tobj, timestep_plus_30s, the_other_one)
# got_intersection_est = False
# if arc_calipso.intersects(arc_the_other_one):
# got_intersection_est = True
got_intersection_acurate = False
arc_calipso_vector = get_arc_vector(tobj, timestep_plus_30s, calipso, arc_len_min)
arc_the_other_one_vector = get_arc_vector(tobj, timestep_plus_30s, the_other_one, arc_len_min)
# Approximate tracs with one arc each arc_len_min minutes.
# For each pair of arcs check if they intersect.
# There is atmost one intersection quit when we find it.
for arc_calipso in arc_calipso_vector:
for arc_the_other_one in arc_the_other_one_vector:
if arc_calipso.intersects(arc_the_other_one):
got_intersection_acurate = True
if got_intersection_acurate:
break
if got_intersection_acurate:
break
if got_intersection_acurate:
intersect = arc_calipso.intersection(arc_the_other_one)
point = (math.degrees(intersect.lon),
math.degrees(intersect.lat))
nextp = the_other_one.get_next_passes(tobj - timedelta(seconds = 60*60),
# SNO around tobj check for passes between +- one hour.
# So that wanted pass for sure is next pass!
2, # Number of hours to find overpasses
point[0],
point[1],
0)
if len(nextp) > 0:
riset, fallt, maxt = nextp[0]
else:
print("No next passes found for, probably a bug!")
tobj = tobj + dtime
continue
nextp = calipso.get_next_passes(tobj - timedelta(seconds = 60*60),
2,
point[0],
point[1],
0)
if len(nextp) > 0:
riset, fallt, maxt_calipso = nextp[0]
else:
print("No next passes found for, probably a bug!")
tobj = tobj + dtime
continue
# Get observer look from Norrkoping to the satellite when it is
# in zenith over the SNO point:
azi, elev = the_other_one.get_observer_look(maxt,
NRK_LON,
NRK_LAT,
NRK_ALT)
isNorrk = (elev > 0.0)
calipsosec = (int(maxt_calipso.strftime("%S")) +
int(maxt_calipso.strftime("%f"))/1000000.)
sec = (int(maxt.strftime("%S")) +
int(maxt.strftime("%f"))/1000000.)
tdelta = (maxt_calipso - maxt)
tdmin = (tdelta.seconds + tdelta.days * 24*3600) / 60.
if abs(tdmin) < minthr:
print(" " +
str(maxt_calipso.strftime("%Y%m%d %H:%M")) +
"%5.1fs" % calipsosec + " "*5 +
str(maxt.strftime("%Y%m%d %H:%M")) +
"%5.1fs" % sec +
" "*6 + "(%7.2f, %7.2f)" % (point[1], point[0]) +
" " + "%4.1f min" % (tdmin) + " " + str(isNorrk)
)
tobj = tobj + timestep_double
if tobj - tobj_tmp > timedelta(days=1):
tobj_tmp = tobj
print(tobj_tmp.strftime("%Y-%m-%d"))
LOG.debug(tobj_tmp.strftime("%Y-%m-%d"))
