geolib.py 文件源码

def wgs84_to_egm96(dem_ds, geoid_dir=None):
    from pygeotools.lib import warplib

    #Check input dem_ds - make sure WGS84

    geoid_dir = os.getenv('ASP_DATA')
    if geoid_dir is None:
        print("No geoid directory available")
        print("Set ASP_DATA or specify")

    egm96_fn = geoid_dir+'/geoids-1.1/egm96-5.tif' 
    try:
        open(egm96_fn)
    except IOError:
        sys.exit("Unable to find "+egm96_fn)
    egm96_ds = gdal.Open(egm96_fn)

    #Warp egm96 to match the input ma
    ds_list = warplib.memwarp_multi([dem_ds, egm96_ds], res='first', extent='first', t_srs='first') 

    #Try two-step with extent/res in wgs84
    #ds_list = warplib.memwarp_multi([dem_ds, egm96_ds], res='first', extent='intersection', t_srs='last') 
    #ds_list = warplib.memwarp_multi([dem_ds, ds_list[1]], res='first', extent='first', t_srs='first')

    print("Extracting ma from dem and egm96 ds")
    from pygeotools.lib import iolib
    dem = iolib.ds_getma(ds_list[0])
    egm96 = iolib.ds_getma(ds_list[1])

    print("Removing offset")
    dem_egm96 = dem - egm96

    return dem_egm96

#Run ASP dem_geoid adjustment utility
#Note: this is multithreaded