ngageoint / sarpy

A basic Python library to demonstrate reading, writing, display, and simple processing of complex SAR data using the NGA SICD standard.
MIT License
262 stars 87 forks source link

DTED regression for voids with 1.3.59 #569

Open kjurka opened 1 week ago

kjurka commented 1 week ago

This PR changes the handling of voids in DTED: https://github.com/ngageoint/sarpy/pull/535/files#diff-7279fd3655c7d057eae1bb04161914b0bc14ae1decdbc06e706b048021ed9f69L373

I have a DTED dataset that has voids wherever there is water. I can't share that here, but I was able to get non void filled DTED tiles from USGS. https://earthexplorer.usgs.gov/ These have voids not over water, but possibly places where they had processing problems. I pulled n33_w119_3arc_v1.dt1 I set up a correctly laid out dted folder and geoid and then ran the below code:

import numpy as np
from sarpy.io.DEM.DTED import DTEDInterpolator, DTEDList

dted_paths = '/share/home/kjurka@emagsys.com/tmp/void_dem_tiles'
ll = np.array([33.938923,-118.388885])

dted_list = DTEDList(dted_paths)
dem_interpolator = DTEDInterpolator.from_reference_point(ll, dted_list, pad_value=0.05)

print(dem_interpolator.get_elevation_hae(ll[0], ll[1]))

With current 1.3.59, I get an elevation of -32803.00222977538 while with 1.3.58, I get an elevation of -37.00222977537916.

The old code looked bogus to me, replacing voids with -1, but the new code doesn't handle them at all. I'd need to study my real DTED dataset more carefully to see if there are any voids other than for water, but if it was just water, replacing with 0 would be reasonable. Your thoughts or other ways to handle this?

kjurka commented 1 week ago

I misspoke about what the voids in my DTED dataset represent. They are not just water/ocean and appear in large/strange shapes/sizes. So replacing with 0 universally is not a good solution (but, I guess better than -32767 m). In the absence of any other data, I'm not sure what could be done. Perhaps if get_elevation_hae returned NaN or some other marker, I could replace with the SCP height.

pressler-vsc commented 1 week ago

Thanks for pointing this out @kjurka. I agree that fixing the reading of DTED nulls in #535 has some strange interactions with DTEDInterpolator and really uncovers the fact that sarpy.io.DEM.DTED is not currently set up to accommodate them. It's possible null-handling in DTED interpolation will be added as a feature in the future. In the interim, I'd recommend accessing the underlying DTED values (e.g. raw_values = DTEDReader('/path/to/dted')[...]) such that you can then control the interpolation, null-handling, bounds-checking, etc. required for your application and DEM source. Does something like that make sense?

kjurka commented 1 week ago

What we've done for now is eliminate DTED tiles from our dataset that have any voids as it was a relatively small percentage in places that we don't care as much about.

Maybe the missing_error flag should be passed further down to the actual DTED reading. If missing_error is True, fail if there are voids. If missing_error is False, then replaces voids with 0, the same as if the tile wasn't there at all. That would be a little more consistent/sane than the current very negative elevations given.