2

I am carrying out some analysis in python (using Anaconda), I have a binary raster and want to perform a geodesic distance calculation from all areas with a value of 1- I can do this within QGIS using the proximity raster tool. In python I have been using gdal, however the gdal_proximity function that QGIS uses does not seem to be available, I have been using the gdal documentation https://gdal.org/python/ and cannot see it there either.

Is this just unavailable when using gdal in python and, if so, is there an alternative?

0

2 Answers 2

5

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.

1
  • I had to delete the quotes around the "VALUES= " in gdal.ComputeProximity to get the correct output pixel values. Commented Nov 26, 2020 at 1:11
0

I have to do the same operation in one task. I don't think QGIS can do the Geodesic distance if the raster file is using the GCS instead of a projection.

Besides, if a raster is on a global scale, the left and right should meet, and the distance should be optimized there, too.

In order to calculate Geodesic distance, such as great circle distance, the raster file needs to be converted to GCS, and a customized distance function should be used in the compute proximity function.

1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.