search

dbekaert / RAiDER

Raytracing Atmospheric Delay Estimation for RADAR
Apache License 2.0
70 stars 39 forks source link

[BUG] Lack of Availability of HRRR Weather Models Means 3-Date Azimuth Grid Interpolation Fails for GUNW workflow #584

Closed cmarshak closed 11 months ago

cmarshak commented 1 year ago 
INFO:RAiDER:Extent of the input is (xmin, ymin, xmax, ymax): -78.05, 41.60, -74.40, 43.55
Traceback (most recent call last):
  File "/opt/conda/envs/RAiDER/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/opt/conda/envs/RAiDER/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/opt/conda/envs/RAiDER/lib/python3.10/site-packages/RAiDER/cli/__main__.py", line 46, in <module>
    main()
  File "/opt/conda/envs/RAiDER/lib/python3.10/site-packages/RAiDER/cli/__main__.py", line 42, in main
    process_entry_point.load()()
  File "/opt/conda/envs/RAiDER/lib/python3.10/site-packages/RAiDER/cli/raider.py", line 620, in calcDelaysGUNW
    cube_filenames = calcDelays([path_cfg])
  File "/opt/conda/envs/RAiDER/lib/python3.10/site-packages/RAiDER/cli/raider.py", line 388, in calcDelays
    raise NotImplementedError(f'The {interp_method} with {n} retrieved weather model files was not well posed '
NotImplementedError: The azimuth_time_grid with 1 retrieved weather model files was not well posed for the dela current workflow.

is raised from these inputs:

{'execution_started': True,
 'expiration_time': '2023-10-25T00:00:00+00:00',
 'job_id': '00312fde-8f89-4c5c-a0d6-ec2df678f316',
 'job_parameters': {'compute_solid_earth_tide': True,
                    'estimate_ionosphere_delay': True,
                    'frame_id': 16421,
                    'granules': ['S1A_IW_SLC__1SDV_20170210T225922_20170210T225949_015228_018F0A_787A',
                                 'S1A_IW_SLC__1SDV_20170210T225947_20170210T230014_015228_018F0A_10FF'],
                    'secondary_granules': ['S1A_IW_SLC__1SDV_20170129T225922_20170129T225949_015053_018996_BEDD',
                                           'S1A_IW_SLC__1SDV_20170129T225947_20170129T230014_015053_018996_6EA8'],
                    'weather_model': 'HRRR'},
 'job_type': 'INSAR_ISCE_TEST',
 'logs': ['https://hyp3-tibet-jpl-contentbucket-81rn23hp7ppf.s3.us-west-2.amazonaws.com/00312fde-8f89-4c5c-a0d6-ec2df678f316/00312fde-8f89-4c5c-a0d6-ec2df678f316.log'],
 'name': 'track106-n3',
 'priority': 100,
 'processing_times': [8558.771, 399.01],
 'request_time': '2023-08-25T03:38:19+00:00',
 'status_code': 'FAILED',
 'user_id': 'access_cloud_based_insar'}

As usual, here is the input GUNW to the above step:

https://hyp3-tibet-jpl-contentbucket-81rn23hp7ppf.s3.us-west-2.amazonaws.com/00312fde-8f89-4c5c-a0d6-ec2df678f316/S1-GUNW-A-R-106-tops-20170210_20170129-225948-00078W_00042N-PP-cde3-v3_0_0.nc

Hoping this can be quickly modified together.

jlmaurer commented 1 year ago

@cmarshak just also documenting the relevant issue here. get_nearest_wmtimes is the relevant function, and possibly the issue is that only one time gets returned because this GUNW has products that are close to the hour. The key will be to ensure that this case gets properly handled. Some additional few lines similar to what I added in #585 should be sufficient.

bbuzz31 commented 1 year ago

Using the above GUNW, I get past the first bug but see a new one:


/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/s1_azimuth_timing.py:276: RuntimeWarning: invalid value encountered in divide
  wgts_norm = [wgt / wgts_sum for wgt in wgts_masked]
CRITICAL: Weather model contains NaNs!
CRITICAL:RAiDER:Weather model contains NaNs!
Processing slice 1 / 20: -500.0
INFO:RAiDER:Processing slice 1 / 20: -500.0
/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/delay.py:275: RuntimeWarning: invalid value encountered in cast
  nParts = np.ceil(ray_lengths.max((1,2)) / MAX_SEGMENT_LENGTH).astype(int) + 1
Traceback (most recent call last):
  File "/Users/buzzanga/mambaforge/envs/RAiDER/bin/raider.py", line 8, in <module>
    sys.exit(main())
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/cli/__main__.py", line 42, in main
    process_entry_point.load()()
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/cli/raider.py", line 620, in calcDelaysGUNW
    cube_filenames = calcDelays([path_cfg])
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/cli/raider.py", line 393, in calcDelays
    wet_delay, hydro_delay = tropo_delay(
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/delay.py", line 90, in tropo_delay
    ds = _get_delays_on_cube(dt, weather_model_file, wm_proj, aoi, height_levels,
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/delay.py", line 163, in _get_delays_on_cube
    wetDelay, hydroDelay = _build_cube_ray(
  File "/Users/buzzanga/Software_InSAR/RAiDER_git/tools/RAiDER/delay.py", line 279, in _build_cube_ray
    fracs = np.linspace(0., 1., num=nparts)
  File "<__array_function__ internals>", line 200, in linspace
  File "/Users/buzzanga/mambaforge/envs/RAiDER/lib/python3.10/site-packages/numpy/core/function_base.py", line 123, in linspace
    raise ValueError("Number of samples, %s, must be non-negative." % num)
ValueError: Number of samples, -9223372036854775807, must be non-negative.
cmarshak commented 1 year ago

It's clear to me that this occurs when not enough weather files are available. This is not an edge case issue.

Here is the relevant output:

Starting to run the weather model calculation
DEBUG:RAiDER:Starting to run the weather model calculation
Requested date,time: 20161223, 22:19
DEBUG:RAiDER:Requested date,time: 20161223, 22:19
Beginning weather model pre-processing
DEBUG:RAiDER:Beginning weather model pre-processing
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
✅ Found ┊ model=hrrr ┊ product=nat ┊ 2016-Dec-23 22:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
/Users/cmarshak/mambaforge/envs/RAiDER/lib/python3.10/site-packages/herbie/core.py:797: UserWarning: This pattern is interpreted as a regular expression, and has match groups. To actually get the groups, use str.extract.
  logic = df.search_this.str.contains(searchString)
👨🏻‍🏭 Created directory: [/Users/cmarshak/Desktop/control_flow_issue/weather_files/hrrr/20161223]
Note: Returning a list of [12] xarray.Datasets because cfgrib opened with multiple hypercubes.
Number of weather model nodes: 1030617
DEBUG:RAiDER:Number of weather model nodes: 1030617
Shape of weather model: (123, 147, 57)
DEBUG:RAiDER:Shape of weather model: (123, 147, 57)
Bounds of the weather model: 1013693.85/1379693.85/2123479.86/2561479.86 (SNWE)
DEBUG:RAiDER:Bounds of the weather model: 1013693.85/1379693.85/2123479.86/2561479.86 (SNWE)
Weather model: HRRR
DEBUG:RAiDER:Weather model: HRRR
Mean value of the wet refractivity: 12.802516
DEBUG:RAiDER:Mean value of the wet refractivity: 12.802516
Mean value of the hydrostatic refractivity: 188.694321
DEBUG:RAiDER:Mean value of the hydrostatic refractivity: 188.694321

======Weather Model class object=====
Weather model time: 2016-12-23 22:00:00
Latitude resolution: 0.02702702702702703
Longitude resolution: 0.02702702702702703
Native projection: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",262.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]
ZMIN: -100.0
ZMAX: 26000.0
k1 = 0.776
k2 = 0.233
k3 = 3750.0
Humidity type = q
=====================================
Class name: hrrr
Dataset: hrrr
=====================================
A: []
B: []
Number of points in Lon/Lat = 123/147
Total number of grid points (3D): 1030617
=====================================

DEBUG:RAiDER:
======Weather Model class object=====
Weather model time: 2016-12-23 22:00:00
Latitude resolution: 0.02702702702702703
Longitude resolution: 0.02702702702702703
Native projection: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",262.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]
ZMIN: -100.0
ZMAX: 26000.0
k1 = 0.776
k2 = 0.233
k3 = 3750.0
Humidity type = q
=====================================
Class name: hrrr
Dataset: hrrr
=====================================
A: []
B: []
Number of points in Lon/Lat = 123/147
Total number of grid points (3D): 1030617
=====================================

Extent of the weather model is (xmin, ymin, xmax, ymax):-70.38, 43.53, -63.57, 47.88
INFO:RAiDER:Extent of the weather model is (xmin, ymin, xmax, ymax):-70.38, 43.53, -63.57, 47.88
Extent of the input is (xmin, ymin, xmax, ymax): -68.60, 44.75, -64.80, 46.70
INFO:RAiDER:Extent of the input is (xmin, ymin, xmax, ymax): -68.60, 44.75, -64.80, 46.70
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
✅ Found ┊ model=hrrr ┊ product=nat ┊ 2016-Dec-23 23:00 UTC F00 ┊ GRIB2 @ aws ┊ IDX @ aws
/Users/cmarshak/mambaforge/envs/RAiDER/lib/python3.10/site-packages/herbie/core.py:797: UserWarning: This pattern is interpreted as a regular expression, and has match groups. To actually get the groups, use str.extract.
  logic = df.search_this.str.contains(searchString)
Note: Returning a list of [12] xarray.Datasets because cfgrib opened with multiple hypercubes.
Number of weather model nodes: 1030617
DEBUG:RAiDER:Number of weather model nodes: 1030617
Shape of weather model: (123, 147, 57)
DEBUG:RAiDER:Shape of weather model: (123, 147, 57)
Bounds of the weather model: 1013693.85/1379693.85/2123479.86/2561479.86 (SNWE)
DEBUG:RAiDER:Bounds of the weather model: 1013693.85/1379693.85/2123479.86/2561479.86 (SNWE)
Weather model: HRRR
DEBUG:RAiDER:Weather model: HRRR
Mean value of the wet refractivity: 12.805363
DEBUG:RAiDER:Mean value of the wet refractivity: 12.805363
Mean value of the hydrostatic refractivity: 188.734711
DEBUG:RAiDER:Mean value of the hydrostatic refractivity: 188.734711

======Weather Model class object=====
Weather model time: 2016-12-23 23:00:00
Latitude resolution: 0.02702702702702703
Longitude resolution: 0.02702702702702703
Native projection: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",262.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]
ZMIN: -100.0
ZMAX: 26000.0
k1 = 0.776
k2 = 0.233
k3 = 3750.0
Humidity type = q
=====================================
Class name: hrrr
Dataset: hrrr
=====================================
A: []
B: []
Number of points in Lon/Lat = 123/147
Total number of grid points (3D): 1030617
=====================================

DEBUG:RAiDER:
======Weather Model class object=====
Weather model time: 2016-12-23 23:00:00
Latitude resolution: 0.02702702702702703
Longitude resolution: 0.02702702702702703
Native projection: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",262.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",38.5,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]
ZMIN: -100.0
ZMAX: 26000.0
k1 = 0.776
k2 = 0.233
k3 = 3750.0
Humidity type = q
=====================================
Class name: hrrr
Dataset: hrrr
=====================================
A: []
B: []
Number of points in Lon/Lat = 123/147
Total number of grid points (3D): 1030617
=====================================

Extent of the weather model is (xmin, ymin, xmax, ymax):-70.38, 43.53, -63.57, 47.88
INFO:RAiDER:Extent of the weather model is (xmin, ymin, xmax, ymax):-70.38, 43.53, -63.57, 47.88
Extent of the input is (xmin, ymin, xmax, ymax): -68.60, 44.75, -64.80, 46.70
INFO:RAiDER:Extent of the input is (xmin, ymin, xmax, ymax): -68.60, 44.75, -64.80, 46.70
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
Weather model HRRR is available from 2016-07-15 to Present
INFO:RAiDER:Weather model HRRR is available from 2016-07-15 to Present
💔 Did not find ┊ model=hrrr ┊ product=nat ┊ 2016-Dec-23 21:00 UTC F00
ERROR:
No index file was found for None
Download the full file first (with `H.download()`).
You will need to remake the Herbie object (H = `Herbie()`)
or delete this cached property: `del H.index_as_dataframe()`
ERROR:RAiDER:
No index file was found for None
Download the full file first (with `H.download()`).
You will need to remake the Herbie object (H = `Herbie()`)
or delete this cached property: `del H.index_as_dataframe()`
WARNING:
WARNING:RAiDER:

The GUNW I am using is: https://hyp3-tibet-jpl-contentbucket-81rn23hp7ppf.s3.us-west-2.amazonaws.com/fe0ef2be-262b-4031-8e22-99b93df661f8/S1-GUNW-A-R-091-tops-20171230_20161223-221939-00069W_00045N-PP-c9b7-v3_0_0.nc

I can easily just weight whatever weather models are available... however, this could yield unexpected results. Not sure if all the weather model files (and the total used) are recorded. Would be important here. @dbekaert @bbuzz31

At some point it will be helpful to understand the distribution of available HRRR model files that are available over Sentinel-1 catalog. It's producing the bulk of our errors for this step.

There are other options too - we could go through download of say 6 closest dates, download them and then just use the 3 that are closest. This has implications related to how long it takes to download data.

cmarshak commented 1 year ago

You can recreate the error running raider.py ++process calcDelaysGUNW -f S1-GUNW-A-R-091-tops-20171230_20161223-221939-00069W_00045N-PP-c9b7-v3_0_0.nc -m HRRR -interp azimuth_time_grid

cmarshak commented 1 year ago

@jlmaurer and @dbekaert just reiterating this is about weather model availability for particular datetimes. We can solve this easily, but given we are in the midst of processing want to utilize something with long term utility.

Previously, for what is now called center_time interpolation i.e. performing an inverse weighting based on 2 closest dates, using 1 or 2 dates was permitted. For the current workflow, I require 3 dates and if 3 dates are not available then the workflow errors out. This is to ensure there is not unexpected behavior.

Here are some possible solutions to the error reported:

  1. Of three closest date times relative to GUNW acq, perform inverse weighting on what is available (as long as there is at least 1). The case of 2 available dates would be equivalent to the interpolation Brett does with two files (except it would be done on a combination of 2 of 3 closest weather model date times relative to GUNW acq) and the case of 1 is using a given date time within 3 that are closest.
  2. Get 3 closest AVAILABLE weather model date times closest to GUNW acq and perform inverse weighting interpolation on these. We could also restrict our search within a window of time of the GUNW acq to error out if 3 weather model files are NOT available in said window.

There may be other solutions too. Again, not sure what these fixes means in terms of interpolating the tropo correction layers save for the ease of stitching.

jlmaurer commented 1 year ago

@cmarshak @dbekaert I can tell you that option 1 will be the easiest, and option 2 is likely the best IMHO. I'm in Guatemala at the moment but I do have some time this evening, could help quickly implement (1) if desired.

cmarshak commented 1 year ago

@jlmaurer - thank you for offer.

Based on conversation with @dbekaert, @bbuzz31 and @sssangha, all the proposed solutions are not going to be used. 😆

We want to do the following:

Given our extended deadline - there is not the same urgency - our processing deadline is likely going to be extended passed FY. In other words, we will record processing issues and make sure we can remedy them in the next month or so.

cmarshak commented 1 year ago

After this Thursday evening (so beginning Friday), we can also merge PRs that you have postponed thanks to the extension of processing deadline. Thank you for your patience @jlmaurer.

jlmaurer commented 1 year ago

@cmarshak sounds like a good plan!!

cmarshak commented 1 year ago

I wanted to provide small update thanks to @asjohnston-asf and @jhkennedy.

The ingest of GUNWs occurs here: https://github.com/asfadmin/grfn-ingest/blob/test/echo10-construction/src/echo10_construction.py#L134

Currently, we only write to the json here: https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/cli/raider.py#L609 (love the use of setdefault as an aside 😉). There are cases where we are writing nothing as here: https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/cli/raider.py#L590-L596

Therefore, as we move towards adding control flow that does not include tropo layer if HRRR datetimes are not available, we would have to ensure these are properly synced up via tests.

Furthermore, want to higlight a few more pieces:

  1. Because it was not known to Sim and myself regarding HRRR availability, we also request 3 weather model times even if only 2 are needed (only certain edge cases when the S1 acq passes over at a model time is when 3 are needed): https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/cli/raider.py#L261-L268. We may want a few lines to see which times 2 or 3 times are needed.
  2. We want to include metadata in which we report which weather model files are used for provenance/reproducibility
  3. There could potentially be weird edge cases in model times for interpolation and may want to add a buffer here so that if s1 time would need a weather model file that is not available yet, we can also stop: https://github.com/dbekaert/RAiDER/blob/dev/tools/RAiDER/aria/prepFromGUNW.py#L101

Again, should have tests for all 3.

cmarshak commented 11 months ago

Resolved by #597