Here is an example that
- opens vector dataset and loops through each feature,
- gdal.Rasterize each feature
- gdal.ComputeProximity each rasterized feature
You should reproject your vector if they are in geographic lat/lon CRS so the proximity (distance) values are in metres, I used QGIS to reproject this data to UTM Zone 18N.
from osgeo import gdal
gdal.UseExceptions()
out_raster_template = "out_{}.tif"
shape_file = "boroughs_utm.shp"
pixel_size = 10
nodata = -9999
id_field = 'boro_code'
drv = gdal.GetDriverByName("ESRI Shapefile")
shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()
xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()
feat_def = lyr.GetLayerDefn()
lyr.ResetReading()
for feat in lyr:
id = int(feat.GetField(id_field))
tmp_feat = feat.Clone()
out_raster = out_raster_template.format(id)
tmp_fn = '/vsimem/tmp.shp'
tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
tmp_lyr.CreateFeature(tmp_feat)
tmp_feat, tmp_lyr, tmp_ds = None, None, None
out_ds = gdal.Rasterize(out_raster, tmp_fn,
outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
noData=nodata, initValues=nodata,
xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
allTouched=True, burnValues=0)
out_ds = None
src_ds = gdal.OpenEx(out_raster, gdal.OF_RASTER)
dst_ds = gdal.OpenEx(out_raster, gdal.OF_UPDATE)
src_band = src_ds.GetRasterBand(1)
dst_band = dst_ds.GetRasterBand(1)
gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES=0'])
dst_band, src_band, dst_ds, src_ds = None, None, None, None
drv.Delete(tmp_fn)
If you want to rasterize some value from your shapefile to use later, but still do the euc distance to another:
from osgeo import gdal
gdal.UseExceptions()
#Create euclidean distance for each polygon and store "Values"
out_raster_template = "out_{}.tif"
out_proximity_template = "prox_{}.tif"
shape_file = "boroughs_utm.shp"
pixel_size = 10
nodata = -9999
id_field = 'boro_code'
value_field = 'value'
drv = gdal.GetDriverByName("ESRI Shapefile")
shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()
xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()
feat_def = lyr.GetLayerDefn()
lyr.ResetReading()
for feat in lyr:
id = int(feat.GetField(id_field))
val = feat.GetField(value_field)
tmp_feat = feat.Clone()
out_raster = out_raster_template.format(id)
prox_raster = out_proximity_template.format(id)
tmp_fn = '/vsimem/tmp.shp'
tmp_raster = '/vsimem/tmp.tif'
tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
tmp_lyr.CreateFeature(tmp_feat)
tmp_feat, tmp_lyr, tmp_ds = None, None, None
out_ds = gdal.Rasterize(out_raster, tmp_fn,
outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
noData=nodata, initValues=nodata,
xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
allTouched=True, burnValues=val)
out_ds = None
out_ds = gdal.Rasterize(tmp_raster, tmp_fn,
outputType=gdal.GDT_Int32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
noData=nodata, initValues=nodata,
xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
allTouched=True, burnValues=id)
out_ds = None
gdal.Translate(prox_raster, out_raster, creationOptions=["COMPRESS=DEFLATE"])
src_ds = gdal.OpenEx(tmp_raster, gdal.OF_RASTER)
dst_ds = gdal.OpenEx(prox_raster, gdal.OF_UPDATE)
src_band = src_ds.GetRasterBand(1)
dst_band = dst_ds.GetRasterBand(1)
gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES={id}'])
dst_band, src_band, dst_ds, src_ds = None, None, None, None
drv.Delete(tmp_fn)
Best Answer
For those battling the same issue in the future: I came up with a new solution, that I posted on GitHub.
TL;DR
: This is the codehttps://gist.github.com/MarByteBeep/e688759bc051a59a5651e0ecd45240ff
The long read:
Based on work done by Spatial Thoughts I created a simple algorithm that iteratively created new features and stores the distance to the original feature in a
BURN
field. These distances are calculated using Azimuthal Equidistant projection. These new features can then be used to create a more accurate proximity map.Let me provide an example. Consider the following features (military objects around some local town)
Next I run the proximity script, with the following settings
Buffer Distance
specifies the distance the buffer is drawn around a feature. In this case I defined a distance of 5000 meters.Segments
specifies how well curves are approximated. Lower is more jaggy.Iterations
specifies how many iso-distances are generated.Rasterize render resolution
specifies the rasterizer resolution in degrees.You can deselect
Open output file after running algorithm
of theProximity Vector layer
, as that layer is mostly an intermediate layer. But for the sake of explaining what's going on I'll leave it switched on.When hitting run two layers,
Proximity Vector layer
andRasterized
, are generated.First let's focus on the
Proximity Vector layer
:When inspecting one of the structures, you can see that all rings contain different
BURN
values, each representing the distance to the original feature. Obviously I could also have named the field:DISTANCE
. But whatever :)The other layer
Rasterize
is the one we're after:When sampling this layer, the value in the band will represent the distance in meters.
Code available below as a GitHub gist. Any optimization suggestions or other ideas are more than welcome!
https://gist.github.com/MarByteBeep/e688759bc051a59a5651e0ecd45240ff
After running 'the long road' example posted in the original question with the new algorithm, I ended up with something more intuitively correct: