CRS hard coded in Geoid conversion

[brentwilder]

Issue:

• CRS is hard coded in the geoid conversion (for the reference, snow-free DEM).
• And also, Ames Stereo Pipeline was not automatically saving the ellipsoid CRS, causing the need for knowing the ellipsoid equivalent CRS after the transformation.

Solution:

Honestly, we should just move this whole piece out of the main pipeline and just keep it as a helper / utility function. If someone needs to transform their vertical datum of the reference DEM to match their lidar, we can just have a separate python file to do this.

Even for our use case, we have not used this once! And it's because we have been using the same ellipsoid DEM for all runs.

    if dem_is_geoid is True:
        # ASP needs NAVD88 conversion to be in NAD83 (not WGS84)
        nad83_dem = join(result_dir, 'demNAD_tmp.tif')
        gdal_func = join(asp_dir, 'gdalwarp')
        cl_call(f'{gdal_func} -t_srs EPSG:26911 {in_dem} {nad83_dem}', log)
        # Use ASP to convert from geoid to ellipsoid
        ellisoid_dem = join(result_dir, 'dem_wgs')
        geoid_func = join(asp_dir, 'dem_geoid')
        cl_call(f'{geoid_func} --nodata_value -9999 {nad83_dem} \
                --geoid NAVD88 --reverse-adjustment -o {ellisoid_dem}', log)
        # Set it back to WGS84
        ref_dem = join(result_dir, 'ellipsoid_DEM.tif')
        cl_call(f'{gdal_func} -t_srs EPSG:32611 {ellisoid_dem}-adj.tif {ref_dem}', log)

        # check for success
        if not exists(ref_dem):
            raise Exception('Conversion to ellipsoid failed')

        log.info('Merged DEM converted to ellipsoid per user input')

    else:
        # cl_call('cp '+ in_dem +' '+ ref_dem, log)
        ref_dem = in_dem
        log.info('Merged DEM was kept in original ellipsoid form...')

[shadoneel]

Brent! Thanks for looking at this. THis use case IS really important, and we should figure out how to have a clear way to use the geoid conversion. We are expanding use into several new domains, and when establishing a new use case, youi dont often have a good reference DEM. FOr example, we have one we are establishing that USGS data is available, but needs cleaning and flight strip mosaic'ing. Ice roads is perfect, except the USGS data is geoid referenced, and we cant do a DEM filter or co-registration because there is no other DEM. SRTM is too coarse given the steep terrain, it erases large regions. Ideally, we would have some good guidance on a geoid-ellipsoid transform, some output that it worked, and then flags to disable DEM filter and co-registration and differencing. In concept, i think these would be super easy to modify the code to do (just some catch statements), but those (like me) who have too clean of hands would spend enormous amounts of time trying to figure it out. Make sense? Can you see a path to getting the functionality for these kinds of operations? Thanks Shad

…

On Wed, Jan 15, 2025 at 2:53 PM Brent Wilder ***@***.***> wrote: Issue: - CRS is hard coded in the geoid conversion. - And also, Ames Stereo Pipeline was not automatically saving the ellipsoid CRS, causing the need for knowing the ellipsoid equivalent CRS after the transformation. Solution: Honestly, we should just move this whole piece out of the main pipeline and just keep it as a helper / utility function. If someone needs to transform their vertical datum of the reference DEM to match their lidar, we can just have a separate python file to do this. Even for our use case, we have not used this once! And it's because we have been using the same ellipsoid DEM for all runs. if dem_is_geoid is True: # ASP needs NAVD88 conversion to be in NAD83 (not WGS84) nad83_dem = join(result_dir, 'demNAD_tmp.tif') gdal_func = join(asp_dir, 'gdalwarp') cl_call(f'{gdal_func} -t_srs EPSG:26911 {in_dem} {nad83_dem}', log) # Use ASP to convert from geoid to ellipsoid ellisoid_dem = join(result_dir, 'dem_wgs') geoid_func = join(asp_dir, 'dem_geoid') cl_call(f'{geoid_func} --nodata_value -9999 {nad83_dem} \ --geoid NAVD88 --reverse-adjustment -o {ellisoid_dem}', log) # Set it back to WGS84 ref_dem = join(result_dir, 'ellipsoid_DEM.tif') cl_call(f'{gdal_func} -t_srs EPSG:32611 {ellisoid_dem}-adj.tif {ref_dem}', log) # check for success if not exists(ref_dem): raise Exception('Conversion to ellipsoid failed') log.info('Merged DEM converted to ellipsoid per user input') else: # cl_call('cp '+ in_dem +' '+ ref_dem, log) ref_dem = in_dem log.info('Merged DEM was kept in original ellipsoid form...') — Reply to this email directly, view it on GitHub <#19>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA22DJNXRYG2T4RAEZ6HVL32K3KFLAVCNFSM6AAAAABVIGYBVOVHI2DSMVQWIX3LMV43ASLTON2WKOZSG44TCMBRGE2TINQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[brentwilder]

@shadoneel yeah absolutely we can keep this -g flag and just improve it to allow for other CRS. I can try working on that this/next week!