Replies: 8 comments 7 replies
-
|
try this (I actually almost had this ready from a GPS class I took): from datetime import datetime
import numpy as np
import pyproj
from lxml import etree as ET
from dateutil.parser import parse
ECEF_TO_LL_TRANSFORMER = pyproj.Transformer.from_crs(
{"proj": "geocent", "ellps": "WGS84", "datum": "WGS84"},
{"proj": "latlong", "ellps": "WGS84", "datum": "WGS84"},
)
def ecef2lla(xyz, degrees=True):
"""Transform a set of [X, Y, Z] ecef points into lat, lon, alt
Inputs:
xyz (ndarray): (N, 3) array of points. each row is [X, Y, Z]
Units: [degrees, degrees, meters] if degrees is True, or
[radians, radians, meters] if degrees is False
Returns: ndarray ()
"""
if xyz.ndim < 2:
xyz = xyz.reshape(1, -1)
if xyz.shape[1] != 3:
raise ValueError("xyz must be (N, 3)")
radians = not degrees
lon_arr, lat_arr, alt_arr = ECEF_TO_LL_TRANSFORMER.transform(
*(xyz.T), radians=radians
)
return np.array([lat_arr, lon_arr, alt_arr]).T
def get_nearest_osv(image_dt, orbit_file):
"""Download/shrink a precise orbit file to just the elements."""
# <List_of_OSVs count="9361">
# <OSV>
# <TAI>TAI=2018-09-30T23:00:19.000000</TAI>
# <UTC>UTC=2018-09-30T22:59:42.000000</UTC>
# <UT1>UT1=2018-09-30T22:59:42.048478</UT1>
# <Absolute_Orbit>+23934</Absolute_Orbit>
# <X unit="m">-1963417.561738</X>
# <Y unit="m">4306539.928523</Y>
# <Z unit="m">5250491.230731</Z>
# <VX unit="m/s">-496.824972</VX>
# <VY unit="m/s">5770.263110</VY>
# <VZ unit="m/s">-4906.500176</VZ>
# <Quality>NOMINAL</Quality>
# </OSV>
tree = ET.parse(orbit_file)
osv_list_elem = tree.find("Data_Block/List_of_OSVs")
for osv in osv_list_elem:
t_orbit = datetime.fromisoformat(osv[1].text[4:])
if t_orbit > image_dt:
x = float(osv[4].text)
y = float(osv[5].text)
z = float(osv[6].text)
return np.array([x, y, z])
raise ValueError("No orbit found for this image")
# Set up simple CLI to take datetime and orbit file
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("datetime", help="Image datetime")
parser.add_argument("orbit_file", help="Precise orbit file")
args = parser.parse_args()
image_dt = parse(args.datetime)
orbit_file = args.orbit_file
xyz = get_nearest_osv(image_dt, orbit_file)
lla = ecef2lla(xyz)
print("Lat,Lon,Alt")
for row in lla:
print(",".join(map(str, np.around(row, 3))))on a test file: so it said the satellite was at Checking on the S1B orbit file, we get something about 180 degrees away (opposite side of earth) |
Beta Was this translation helpful? Give feedback.
-
|
Most robust way of determining which orbit to use is to look at sensor name |
Beta Was this translation helpful? Give feedback.
-
|
Piyush, the metadata is wrong in the GUNW 😢 - not by choice of course - so if you need all the SLCs for that generated the igram, I don't think this would work particularly with older GUNW versions. That said, @bbuzz31 - I think there is one reference/secondary scene that is correct within a GUNW consistently (you only need one scene to get the correct orbit file I think) - I asked a similar question here: #391. This should be an attribute within the GUNW dataset. Scott's way of geolocating is nice. I like this. I could potentially see some edge cases when the two sensors cross the same area in a given day (sensing time is not readily available in the GUNW either). But maybe I am wrong about this. I think you need exact sensing time to get this to work - you would have to decouple ascending/descending tracks from the orbit file (fortunately, this is in the GUNW file name). If you are not opposed to hitting CMR, you can always get the exact SLCs of a GUNW (these are always correct) using Andrew Johnston's functions that I documented here: https://github.com/ACCESS-Cloud-Based-InSAR/access-assorted-notebooks/blob/main/marshak/2022_05/0_GUNW_ID_to_Input_Scenes.ipynb (using ASF-search won't work because it only stores metadata for the primary scene - yea that's annoying too). Then you can use Piyush's suggestion, Scott's downloader, or s1_reader for downloading the orbits. |
Beta Was this translation helpful? Give feedback.
-
|
Beta Was this translation helpful? Give feedback.
-
|
Just one name for ref / SLC should work with the 2 min window that we are using by default |
Beta Was this translation helpful? Give feedback.
-
|
@cmarshak Can you remind me what was the bug of the incorrect meta-data of labeling of the reference and secondary SLC names in the project. You mentioned before that at least 1 was always correct, can you elaborate. Can you also confirm its the case for all GUNW product versions? @piyushrpt. My fallback lacking the meta-data information on the orbit would be; fetch both orbits for S1A and S1B; then based on centroid of the coordinate which we have of the product (simple gdal stats) do a geo2rdr using ISCE3 module and then you would get an error return for the one that the location does not match. We could implement that trajectory for all the GUNW product versions that have it potentially incorrect. |
Beta Was this translation helpful? Give feedback.
-
|
If you have the SLC granule name, you have both the time and sensor. That should be sufficient. |
Beta Was this translation helpful? Give feedback.
-
Beta Was this translation helpful? Give feedback.
-
|
A single missing SLC should not be an issue here at all - this is only about 30-40 secs in time. If we use a window of few minutes like we currently are having a nominal granule name should be more than sufficient |
Beta Was this translation helpful? Give feedback.
-
|
Confirmed now with Charlie that there is always a correct SLC listed, but we would just be missing a complete set of SLC. Agreed with @piyushrpt that that is sufficient to now fetch the orbit directly as it has enough information and you have sensor name and time already in it. |
Beta Was this translation helpful? Give feedback.
-
|
There aer couple of things you should consider / test regarding InSAR:
|
Beta Was this translation helpful? Give feedback.
-
|
Very interesting... learned something new. Curious how this will be handled in raider and how it is done elsewhere. I guess for generating an interferogram, you can't do this (or the weird geometry is taken care of in co-registration or something like that). Ok - will be quiet now. |
Beta Was this translation helpful? Give feedback.
-
|
We would be using single orbit and calculate the delays for the reference and secondary from it (i.e. before after and then interferometric). So should indeed cancel out. Tagging @bbuzz31 for his awareness. |
Beta Was this translation helpful? Give feedback.
-
|
Looks like we are all good on this discussion so it can be closed. |
Beta Was this translation helpful? Give feedback.
-
|
Thanks for all the input everyone! |
Beta Was this translation helpful? Give feedback.
-
Background
To support #434 it's necessary to grab the correct orbit file for the dates in the GUNW. Some of the GUNW products have incorrect SLCs stored in the metadata, for now precluding use of https://github.com/opera-adt/s1-reader.
We're currently using the https://github.com/scottstanie/sentineleof which allows us to download orbit files based on date and time. Because the orbit files are two delays long we can potentially get S1A and S1B orbit files (e.g.,
eof -d 2018-10-01downloads S1B_OPER_AUX_POEORB_OPOD_20210314T200444_V20180930T225942_20181002T005942.EOF and S1A_OPER_AUX_POEORB_OPOD_20210309T100229_V20180930T225942_20181002T005942.EOF.Since we have the GUNW, we have geo information which we should be able to use to figure out the correct orbit file.
Any ideas on how to use the geo lat/lon information in the GUNW along with RAiDER tools to figure out which is the correct orbit file?
Hoping @piyushrpt might have some ideas.
Also may be of interest to @scottstanie @jlmaurer @dbekaert
Beta Was this translation helpful? Give feedback.
All reactions