pytroll / pyorbital

Orbital and astronomy computations in python
http://pyorbital.readthedocs.org/
GNU General Public License v3.0
224 stars 77 forks source link

get_lonlatalt doesn't give correct results for Terra / Aqua #70

Open espg opened 3 years ago

espg commented 3 years ago
from datetime import datetime, timezone
from pyorbital import orbital
from pyhdf.SD import SD, SDC

t = datetime(2005, 6, 16, 1, 20, 0, 0, tzinfo=timezone.utc)
terra = orbital.Orbital("EOS-TERRA")

terra.get_lonlatalt(t)

Above code yields:

(166.84878227675836, 60.28586375275336, 718.455146236583)

Terra MODIS is a nadir looking instrument, so those coordinates (or relatively close to them) should be present in the data from the sensor... but...

pass1 = SD('MOD07_L2.A2005167.0115.061.2017291020316.hdf', SDC.READ)
pass2 = SD('MOD07_L2.A2005167.0120.061.2017291020310.hdf', SDC.READ)
lats1 = pass1.select('Latitude')
lons1 = pass1.select('Longitude')
lats2 = pass2.select('Latitude')
lons2 = pass2.select('Longitude')

print(np.max(lats1[:]),np.min(lats1[:]),np.max(lats2[:]),np.min(lats2[:]))

returns:

(-54.213318 -77.61293 -38.029873 -59.76777)

All the latitudes are in the southern hemisphere, while the pyorbital thinks that the satellite is close to the arctic circle in the northern hemisphere, so that's not great...

mraspaud commented 3 years ago

Hi @espg , thanks for reporting this issue.

Calling Orbital without providing specific TLEs does go and fetch the current TLEs from celestrak, meaning the most recent ones. TLEs are supposed to be valid for max a week, since the SGP4 model is a simplification of reality. So trying to extrapolate the position of the satellite over 15 years is stretching it a bit ;) No wonder it's all wrong...

Celestrak.com is providing some historical TLEs, but I think they stop in 2004. However, if you contact them, they can probably send you the rest of the archive.

espg commented 3 years ago

I requested/downloaded the TLE's for both satellites from launch till the beginning of this month...I'm wondering if it makes sense to open a PR that checks for if the date is earlier than some cutoff (a month? 2020?), and use the TLE archive if that's the case, but download current from celestrak if the time is more recent. The text files aren't that large for the archive (1.7 meg's for Aqua to get all TLE's from launch to now). @mraspaud thoughts?

pnuu commented 3 years ago

I thought we already did select the closest TLE in time, but apparently not: https://github.com/pytroll/pyorbital/blob/master/pyorbital/tlefile.py#L197-L201

mraspaud commented 3 years ago

@espg good idea! As a first step, we could issue a loud warning when TLEs are too old indeed. As for getting the data from the archive, the problem is that the archive isn't complete, so it might be difficult, but PRs are definitely very welcome!

espg commented 3 years ago

So... I do have a fix for this, but right now it's fairly klugey. @mraspaud you're right in the public archive being incomplete, but I have requested the historic data, and have a full archive from launch to Jan. 1, 2020. I'm having some interface issues; initially I created single text files containing all the TLE's, but calling pyorbital with a TLE file argument expects just a single TLE-- if there are multiple TLE's, it will only read the first two lines. For a workaround, I took the full archive for both TERRA and AQUA, and split the TLE single text file into per-day TLE's from the archive. The celestrek archive has 3 to 4 TLE's per day, so each of these text files isn't a single TLE, but that doesn't matter, because reading the first TLE for a given day is plenty accurate, and doesn't cause the pyorbital interface to break. In short, right now there is a translator wrapper that parses the datetime object before creating the satellite query object, and translates what the correct TLE textfile is based of a year/day of year naming scheme.

This isn't great. Ideally, when creating a satellite object to query, multiple TLE's could be specified (i.e., an archive of them), and at run time get_lonlatalt(time) would use the time argument to decide which TLE inside of the text file to parse. If the time argument is outside of the archive (newer than the newest/last TLE in the archive), then just grab the current TLE like is currently done. The parsing part isn't bad, since there's already predictable time/date info in the TLE itself to find that is position based. The archives are pretty small-- 3.3 MB for 20+ years of 3-4 TLE's a day on the TERRA sensor, as a text file...compressed as a tar.gz is under a megabyte. If data size is a worry, the archive could be thinned to a single TLE per day, or even 1 per week-- and if there's interest, pyorbital could just host this archive internally in the codebase so it's installed with the library, and select the closest TLE the way that @pnuu and myself expected was already happening. Less ideally, if there isn't interest in pyorbital carrying the archive, the same logic for parsing multiple TLE's in a text file would still enable users to just specify a single archive file... the celestrek standard data request form produces just a simple text file of concatenated TLE's, so that format seems standard for what a user might expect as input. Another option (which is not my preferred solution) is if there's a way to make orbital.Orbital accept a text string instead of a filename or url; that would at least enable external archives on the user side to only have a single file for the archive to track and search/parse at runtime, instead of thousands of small text files or having to create temporary throw-away files per query.

I'm happy to put some sweat and work into making the first and/or second solutions work and submit a PR to do this... but since the current design pattern is to fix the TLE at init of the orbital object, changing things to parse at runtime based on utc time will touch a lot of other parts of the object, so I may need some feedback and support as I work out the implementation. Thoughts?

mraspaud commented 3 years ago

@espg wow, thanks for looking into this!

About Orbital not looking through an archive file, this comes from the fact that it's original purpose was to work with individual TLE files, but using archives is of course a natural enhancement with I support totally. For the record, we do have some code in pygac for looking at the nearest tle lines matching a date here: https://github.com/pytroll/pygac/blob/master/pygac/reader.py#L659 But this is clearly not in the right place. So, I think your proposal of having the tle lookup done at runtime is good. To keep things clean, maybe the tle file could be wrapped into a small class with a method (eg get_tle_closest_to(time)) that would be called by eg get_lonlatalt when needed?

About having archive files in the pyorbital repo: I'm not really inclined to include the tle archive in the pyorbital repo, for multiple reasons, like:

Again, thanks a lot for putting your time on this, it's highly appreciated!