Geodesic distance from raster in python?

Gdal_proximity is in fact a py file that uses gdal.ComputeProximity(srcband,dstband,options) You can check that .py here: https://www.gdal.org/gdal_proximity.html The source code is here: https://svn.osgeo.org/gdal/trunk/gdal/swig/python/scripts/gdal_proximity.py

As you can see in the last part of the source code it uses gdal.ComputeProximity method. Therefore if you don't want to call Gdal_proximity.py you can use a very simple python example as:

    from osgeo import gdal
    
    src_ds = gdal.Open('DEM.tif')
    srcband=src_ds.GetRasterBand(1)
    dst_filename='output.tif'
    
    drv = gdal.GetDriverByName('GTiff')
    dst_ds = drv.Create( dst_filename,
                         src_ds.RasterXSize, src_ds.RasterYSize, 1,
                         gdal.GetDataTypeByName('Float32'))
    
    dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
    dst_ds.SetProjection( src_ds.GetProjectionRef() )
    
    dstband = dst_ds.GetRasterBand(1)
        
    # In this example I'm using target pixel values of 100 and 300. I'm also using Distance units as GEO but you can change that to PIXELS.
    gdal.ComputeProximity(srcband,dstband,["VALUES='100,300'","DISTUNITS=GEO"])
    srcband = None
    dstband = None
    src_ds = None
    dst_ds = None

You can check documentation links above for more information.