Topics

gdal_proximity


Spencer Gardner
 

Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer


Anderson Roberto da Silva
 

Hello,
First sorry for my bad English.
Maybe I can help you and other people with the same doubt. Currently I' ve been working on a project that I need this functionality.

def proximity(in_raster, out_raster, list_values=[]):
"""Calculate proximity
:param in_raster: File path from input raster
:param out_raster: File path from output raster
:param list_values: List of the values.
If not specified, it's considered 0 value
:return: Return the file path from output raster proximity
"""
if len(list_values) == 0:
values = 'VALUES=' + str(0)
else:
values = 'VALUES=' + ','.join(str(value) for value in list_values)

in_raster_ds = gdal.Open(in_raster, gdal.GA_ReadOnly)
x_size = in_raster_ds.RasterXSize
y_size = in_raster_ds.RasterYSize

raster_band1 = in_raster_ds.GetRasterBand(1)

driver = gdal.GetDriverByName('GTiff')

dest_ds = driver.Create(out_raster, x_size, y_size, 1, gdal.GDT_Float32,
creation_options)

dest_ds.SetProjection(in_raster_ds.GetProjection())
dest_ds.SetGeoTransform(in_raster_ds.GetGeoTransform())

dest_ds_band1 = dest_ds.GetRasterBand(1)

gdal.ComputeProximity(raster_band1, dest_ds_band1, [values,
"DISTUNITS=GEO"])

in_raster_ds, raster_band1, dest_ds = None, None, None

return out_raster

For test:
class TestProximity(unittest.TestCase):
in_raster = path.join(RASTER_FOLDER, 'rodovias_raster.tif')
out_raster = path.join(RASTER_FOLDER, 'rodovias_proximity.tif')

def test_proximity(self):
proximity(self.in_raster, self.out_raster)

def test_proximity_values(self):
proximity(self.in_raster, self.out_raster, [1])

But first, you need to rasterize your roads.
My function to this:
def rasterize_fit_by_layer(geopackage, target_layer, field, grid_size,
out_raster, extend_layer=None):
"""Rasterize vector layer to image only in the layer extents

:param geopackage: Geopackage file
:param target_layer: Name of the layer on geopackage to rasterize
:param field: Name of the field that contains the integer values
:param grid_size: Cell size of the grid raster
:param extend_layer: Name of the layer on geopackage to extents (Optional)
:param out_raster: File path to save output raster
:return: File path from output raster
"""
source = ogr.Open(geopackage)

target_layer_rasterize = source.GetLayer(target_layer)
if target_layer_rasterize is None:
return None

layer_def = target_layer_rasterize.GetLayerDefn()
if layer_def.GetFieldIndex(field) < 0:
return None

spatial_ref = target_layer_rasterize.GetSpatialRef()

if extend_layer is None:
extents = target_layer_rasterize.GetExtent()
else:
layer_extent = source.GetLayer(extend_layer)
if layer_extent is None:
return None
extents = layer_extent.GetExtent()

x_min, x_max, y_min, y_max = extents

cols = int((x_max - x_min) / grid_size)
rows = int((y_max - y_min) / grid_size)

x_res = (x_max - x_min) / float(cols)
y_res = (y_max - y_min) / float(rows)

driver = gdal.GetDriverByName('GTiff')
dest_ds = driver.Create(out_raster, cols, rows, 1, gdal.GDT_Byte,
creation_options)
dest_ds.SetGeoTransform((x_min, x_res, 0, y_max, 0, -y_res))
dest_ds.SetProjection((spatial_ref.ExportToWkt()))

options = ['ATTRIBUTE=' + field]
# options = ['ATTRIBUTE=' + field, 'ALL_TOUCHED=TRUE']
gdal.RasterizeLayer(dest_ds, [1], target_layer_rasterize, None,
options=options)

dest_ds, source = None, None
layer_extent, target_layer_rasterize = None, None
Change for your uses.

Hope this helps.

-- 
Anderson Roberto da Silva
Engenheiro Cartógrafo





On Mon, Jan 4, 2021 at 5:21 PM Spencer Gardner <spencergardner@...> wrote:
Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer


Sean Gillies
 

Hi Spencer,

On Mon, Jan 4, 2021 at 1:21 PM Spencer Gardner <spencergardner@...> wrote:
Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer

You're not overlooking anything. The rasterio project doesn't surface GDALComputeProximity. There are pointers to equivalent functionality in scipy at https://github.com/mapbox/rasterio/issues/531.

Hope this helps,

--
Sean Gillies


Spencer Gardner
 

Thanks. I'll give that a go.


On Wed, Jan 13, 2021 at 7:45 AM Sean Gillies <sean.gillies@...> wrote:
Hi Spencer,

On Mon, Jan 4, 2021 at 1:21 PM Spencer Gardner <spencergardner@...> wrote:
Hi all. I'm using Rasterio for a raster analysis but I need to calculate proximity as is done in gdal_proximity.

Is this implemented somewhere in Rasterio? I don't see anything like it in the docs. Or if not, is there a way to run the GDAL operation on the array? Or even better, I could run a numpy operation on it, but I'm not finding anything in numpy that is suitable.

As an example, I am rasterizing a line feature that represents a segment of road. My analysis needs to account for proximity to the road such that cells touched by the road or immediately adjacent are weighted higher than roads that are 500 ft away. My intent with the proximity analysis was to use results to generate these weights.

Any pointers are much appreciated. Thanks!

Spencer

You're not overlooking anything. The rasterio project doesn't surface GDALComputeProximity. There are pointers to equivalent functionality in scipy at https://github.com/mapbox/rasterio/issues/531.

Hope this helps,

--
Sean Gillies