Add script to get SNOs for two platforms

[adybbroe]

...typically Calipso and an Imager satellite like Metop/NOAA

Signed-off-by: Adam Dybbroe Adam.Dybbroe@smhi.se

• Closes #xxxx
• Tests added
• Tests passed
• Passes flake8 pyorbital
• Fully documented

[coveralls]

Coverage decreased (-7.3%) to 76.609% when pulling 91f9f74 on adybbroe:add-sno-finder-script into 7efe253 on pytroll:master.

[adybbroe]

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 ?

[pnuu]

Few comments, and one change request (the addition of main()). Otherwise looks good to me.

[pnuu]

Encapsulate all this in main() function, and in if __name__ == "__main__": block call that. Can save some time later on when hunting a bug that is caused by polluted namespace.

[adybbroe]

Done!

[pnuu]

How about finding the platforms by satellite ID number? The platform names might not always be in the TLE data, but the ID is.

[adybbroe]

@pnuu Does pyorbital know how to do that, thought it required a name?

[mraspaud]

@adybbroe where are we at with this PR ?

[zxdawn]

@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?

[ninahakansson]

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"))

[adybbroe]

Hi!
@ninahakansson and @zxdawn I am opening a new PR based around Ninas (improved) script. Sorry for this PR being dormant for so long. But I have an actual use case now where I want to find SNOs for MW sounders. So, hopefully, I can get it all the way now. Appreciate any testing and feedback and help once I have the new PR. Should come during the day hopefully!

[adybbroe]

@zxdawn I am doing good progress, but not yet ready to commit. Think that will be done in the morning (Thursday) CET.
If you still have some "truth" I can build some test cases around that? By "truth" I mean some actual data where you have the SNO match in terms of time(s) and lon,lat

[zxdawn]

@adybbroe Nice to hear you're working on that! I don't have data at my hand now. Maybe you can test the CALIPSO and TROPOMI cases above? If it makes sense, then I can try to fetch the data later and find the exact location.

[adybbroe]

So, now there is a new draft PR here #153 taking over from this one. I am closing this one (without merging of course).
I will look at your Tropomi case (almost) next @zxdawn :-)