Cop 30m WGS84 ellipsoidal height correction is failing at the Anti-meridian (Antarctica)

[abradley60]

The bug

WGS84 Ellipsoidal height correction (dst_ellipsoidal_height = True) at the anti-meridian is failing on the antarctic coastline.

To Reproduce

requirements

pip install dem-stitcher
pip install asf-search

Query the coordinates of a scene known to cross the antimeridian in Antarctica using the asf-search package:

import asf_search as asf

scene_list = ['S1A_IW_SLC__1SDV_20191028T131928_20191028T131947_029659_0360CA_A1A0']
results = asf.granule_search(scene_list, asf.ASFSearchOptions(processingLevel='SLC))
geom = results[0].geometry
geom

{'coordinates': [[[176.881088, -78.176201],
   [-177.884048, -77.806282],
   [178.838364, -75.697151],
   [174.300995, -76.01265],
   [176.881088, -78.176201]]],
 'type': 'Polygon'}

Note this scene crosses the antimeridian. However, it is a valid polygon that extends the width of the earth so it is not handled by the dem-stitcher software. If we were to pass the bounds of this to the stitch_dem function it will attempt to download a large number of tiles.

Manually split the geometry to a set of bounds to the left and right of the antimeridian:

bounds_left = (-180, -78.176201, -177.884048, -75.697151)
bounds_right = (174.300995, -78.176201, 180, -75.697151)

Download the left dem using dem-stitcher with dst_ellipsoidal_height = True

from dem_stitcher import stitch_dem

dem_data, dem_meta = stitch_dem(bounds_left,
    dem_name='glo_30',
    dst_ellipsoidal_height=True,
    dst_area_or_point='Point',
    merge_nodata_value=0
    )

with rasterio.open('bounds_left.tif', 'w', **dem_meta) as ds:
    ds.write(dem_data,1)
    ds.update_tags(AREA_OR_POINT='Point')

RuntimeError                              Traceback (most recent call last)
Cell In[134], line 3
      1 from dem_stitcher import stitch_dem
----> 3 dem_data, dem_meta = stitch_dem(bounds_left,
      4     dem_name='glo_30',
      5     dst_ellipsoidal_height=True,
      6     dst_area_or_point='Point',
      7     merge_nodata_value=0
      8     )
     10 with rasterio.open('bounds_left.tif', 'w', **dem_meta) as ds:
     11     ds.write(dem_data,1)

File ~/dem-stitcher/dem_stitcher/stitcher.py:362, in stitch_dem(bounds, dem_name, dst_ellipsoidal_height, dst_area_or_point, dst_resolution, n_threads_reproj, n_threads_downloading, fill_in_glo_30, merge_nodata_value)
    359     zipped_data = list(map(lambda ds: _translate_one_tile_across_dateline(ds, crossing), datasets))
    360     memory_files, datasets = zip(*zipped_data)
--> 362 dem_arr, dem_profile = merge_and_transform_dem_tiles(
    363     datasets,
    364     bounds,
    365     dem_name,
    366     dst_ellipsoidal_height=dst_ellipsoidal_height,
    367     dst_area_or_point=dst_area_or_point,
    368     dst_resolution=dst_resolution,
    369     num_threads_reproj=n_threads_reproj,
    370     merge_nodata_value=merge_nodata_value,
    371     n_threads_for_reading_tile_data=n_threads_downloading,
    372 )
    374 # Close datasets
    375 list(map(lambda dataset: dataset.close(), datasets))

File ~/dem-stitcher/dem_stitcher/stitcher.py:198, in merge_and_transform_dem_tiles(datasets, bounds, dem_name, dst_ellipsoidal_height, dst_area_or_point, dst_resolution, num_threads_reproj, merge_nodata_value, n_threads_for_reading_tile_data)
    196 if dst_ellipsoidal_height:
    197     geoid_name = DEM2GEOID[dem_name]
--> 198     dem_arr = remove_geoid(
    199         dem_arr,
    200         dem_profile,
    201         geoid_name,
    202         dem_area_or_point=dst_area_or_point,
    203     )
    205 if dst_resolution is not None:
    206     dem_profile_res = update_profile_resolution(dem_profile, dst_resolution)

File ~/dem-stitcher/dem_stitcher/geoid.py:72, in remove_geoid(dem_arr, dem_profile, geoid_name, dem_area_or_point, res_buffer)
     68 assert dem_area_or_point in ['Point', 'Area']
     70 extent = array_bounds(dem_profile['height'], dem_profile['width'], dem_profile['transform'])
---> 72 geoid_arr, geoid_profile = read_geoid(geoid_name, extent=list(extent), res_buffer=res_buffer)
     74 t_dem = dem_profile['transform']
     75 t_geoid = geoid_profile['transform']

File ~/dem-stitcher/dem_stitcher/geoid.py:46, in read_geoid(geoid_name, extent, res_buffer)
     42 extent_l, extent_r = split_extent_across_dateline(extent)
     43 geoid_arr_l, geoid_profile_l = read_raster_from_window(
     44     geoid_path, extent_l, extent_crs, res_buffer=res_buffer
     45 )
---> 46 geoid_arr_r, geoid_profile_r = read_raster_from_window(
     47     geoid_path, extent_r, extent_crs, res_buffer=res_buffer
     48 )
     49 res_x = geoid_profile_l['transform'].a
     50 if crossing == 180:

File ~/dem-stitcher/dem_stitcher/rio_window.py:161, in read_raster_from_window(raster_path, window_extent, window_crs, res_buffer)
    158 with rasterio.open(raster_path) as ds:
    159     src_profile = ds.profile
--> 161 window = get_window_from_extent(src_profile, window_extent, window_crs, res_buffer=res_buffer)
    163 with rasterio.open(raster_path) as ds:
    164     arr_window = ds.read(window=window)

File ~/dem-stitcher/dem_stitcher/rio_window.py:105, in get_window_from_extent(src_profile, window_extent, window_crs, res_buffer)
     99         warn(
    100             f'Requesting extent beyond raster bounds of {list(src_bounds)}'
    101             f'. Shrinking bounds in raster crs to {window_extent_r}.',
    102             category=RuntimeWarning,
    103         )
    104 else:
--> 105     raise RuntimeError('The extent you specified does not overlap' ' the specified raster as a Polygon.')
    107 corner_ul, corner_br = get_indices_from_extent(
    108     src_profile['transform'], window_extent_r, shape=src_shape, res_buffer=res_buffer
    109 )
    110 row_start, col_start = corner_ul

RuntimeError: The extent you specified does not overlap the specified raster as a Polygon.

Note, if we set dst_ellipsoidal_height = False the dem will download as expected.

In the case of the right side of the antimeridian with dst_ellipsoidal_height = True , the function will work but the dem has some missing data near the antimeridian compared to the case where dst_ellipsoidal_height = False (see image below)

dem_data, dem_meta = stitch_dem(bounds_right,
    dem_name='glo_30',
    dst_ellipsoidal_height=True,
    dst_area_or_point='Point',
    merge_nodata_value=0
    )

with rasterio.open('bounds_right.tif', 'w', **dem_meta) as ds:
    ds.write(dem_data,1)
    ds.update_tags(AREA_OR_POINT='Point')

I'm guessing this has something to do with the coverage of the geoid?

Additional context

I am processing sentinel-1 scenes over Antarctica and acquiring the cop 30m using the dem-stitcher software. This issue has occurred at the antimeridian, elsewhere the dem acquisition is working as expected

[abradley60]

I'm at the end of my workday so just noting my findings here. There appears to be two separate issues:

1. The nodata on the right side of the AM:

I think the missing data on the right side is due to the extent of the geoid that is being applied. The s3://aria-geoid/egm08_25.tif (referenced in the dem-stitcher code here) has the following extent:

The EGM2008 2.5' geoid model that can be accessed here appears to have identical values but a larger extent with an extra pixel at the border of the right side of the antimeridian (180)

Therefore, I think the values in the dem to the right of the antimeridian are being coverted to nodata, as there is missing data in the egm08_25.tif geoid file. Perhaps the geoid file should be updated?

2. Failing on the left side of the antimeridian:

I think this is due to the antimeridian crossing check here. Because the dem extent is slightly outside of the WGS84 range (+180,-180), it attempts to split the extent across the antimeridian, resulting in two extents, one of which is very small. Given the extent of the geoid is (-180.0208333333333428, 180.0208333333332860) these values could be used here instead when the function is called here.

i.e. for the geoid dateline crossing check, change:

if (xmin > -180) and (xmax < 180):

to

if (xmin > -180.0208333333333428) and (xmax < 180.020833333333286):

This solved the issue for me with the bounds_left = (-180, -78.176201, -177.884048, -75.697151)

Let me know if this is not clear, I can raise a PR tomorrow with some code changes that might be easier to follow. Thanks!

[cmarshak]

Thank you so much for the detailed explanation as well as the reproducible calls. Will have to look at it over the course of this week.

Sounds like Issue #100 might solve this if the buffered geoid does resolve this in this anti-meridian area.

At some point today - I am going to upload both the 2.5 and 1 degree geoids from agi-soft to our s3 bucket and share the links here.

The quickest fix would be to swap out the geoid hardcoded here: https://github.com/ACCESS-Cloud-Based-InSAR/dem-stitcher/blob/dev/dem_stitcher/geoid.py#L14-L19

I still think there might be a bug in the code though with respect to this based on what you are saying.

This probably has all come to pass because of the issue here: #96

So to be clear, we were using the 1 degree geoid from agisoft, but it went down previously and that's "no bueno".

Thank you again.

[cmarshak]

It's taking some time on my part - sorry for the delay.

I uploaded the geoids (both 1 and 2.5 degrees) to the ARIA bucket.

I will create a PR and inspect it tomorrow (hopefully) or maybe Thursday.

They file names are:

us_nga_egm2008_25_4326__agisoft.tif
us_nga_egm2008_1_4326__agisoft.tif

Where I used rio warp <input> <output> --dst-crs EPSG:4326 as the software expects the geoid in 4326 (the bounds did not change).

Feel free to check too.

Thanks again.

[abradley60]

Hi @cmarshak no problem at all, I have a working implementation in the meantime.

There are a number of other DEM handling steps our team are doing in order to get the desired format for our sentinel-1 processing. I am thinking it could make sense to add these directly to the dem-stitcher software. Other users may have the same requirements. I might raise a feature suggestion over the next week or two to discuss if your team is open to this?

[cmarshak]

That's good to hear. Thanks for your patience.

We use this software for our ARIA-S1-GUNWs, which are actively being processed and distributed via the DAAC. We don't have to worry about the anti-meridian issue you raise because ISCE2, which is what we use for the INSAR processing, does not work on the antimeridian anyways :) - see isce-framework/isce2#634

Definitely want to address exposing the geoid parameter as mentioned in the issue above to allow other users to BYOG ('bring your own geoid').

The other bugs will take some time to address.

[cmarshak]

I had some time to fix some of the outstanding issues after months of neglect. I (think) I now understand the issue.

It's because the read from both "ends" of the dateline causes np.nan to corrupt the data if that data is unavailable.

For concreteness:

with rasterio.open('https://aria-geoid.s3.us-west-2.amazonaws.com/egm08_25.tif') as ds:
    b = ds.bounds
# b =  BoundingBox(left=-180.02083333333334, bottom=-90.02083333333334, right=179.97916666666666, top=90.02083333333331)

with rasterio.open('https://aria-geoid.s3.us-west-2.amazonaws.com/us_nga_egm2008_1_4326__agisoft.tif') as ds:
    b_agisoft = ds.bounds
#b_agisoft = BoundingBox(left=-180.00833333333333, bottom=-90.00833333333331, right=180.00833333333333, top=90.00833333333334)

When I took care of the dateline initially, I was using the agi-soft geoid (2nd url) which has >=1 pixel buffer around both sides of the dateline. However the AGI soft geoid went down for an extended period (see #96) and swapped it out with the first. The AGI soft has since gone back up and it is now mirrored on our ARIA bucket.

Indeed if there is missing data near the dateline, the process for dealing with dateline crossing (reading from both ends of the dateline) causes issues because all nodata is cast to np.nan, which when added to floats remains np.nan (this is how geoid is "removed"). This results in your np.nan.

So, the way I am planning on dealing with this is:

1. Using the agisoft EGM 2008 geoid for GLO-30
2. If a different geoid is supplied to raise a Userwarning if the data doesn't extend past the dateline (recommend 1 pixel buffer).

Will likely release these fixes in the coming weeks. Thanks again for your thorough documentation of this issue!