# Buffering around raster using gdal and numpy?

I have a stream raster with stream cells containing stream-order values. I need a buffer of 1000m around the streams and it would be good if the values in the buffer raster contain the stream-order to which it is close to.

I am using gdal and numpy modules to handle raster.

How can I do this?

• A good question should include a degree of research and attempt. What have you tried already? What happens when you try it? What have you turned up in your research, and how does that not work for you? Please edit your question to include more detail.
– Midavalo
Commented Jul 31, 2017 at 20:02
• Why the hell you guys down vote? his answer below is awesome! Commented Jul 31, 2017 at 20:36
• @nickves This is a question and answer site. Votes for questions are separate to votes for answers. An excellent answer doesn't mean there was an excellent question.
– Midavalo
Commented Jul 31, 2017 at 20:51
• I thought my initial question was reciprocating adequate information for the answer that I needed. Anyways, I've edited my question a bit and hope it looks better now. I was really surprised to see the downvote because when I search around the internet, there was no forum giving me distance buffering algorithm and I thought this would be a good question. Commented Aug 1, 2017 at 4:17

I was able to build my own algorithm for this and it's working like a charm.

Since I didn't find this anywhere when I googled it, I'm posting my code here in case someone needs it.

from osgeo import gdal
import numpy as np,sys

def raster_buffer(raster_filepath, dist=1000):
"""This function creates a distance buffer around the given raster file with non-zero values.
The value in output raster will have value of the cell to which it is close to."""
d=gdal.Open(raster_filepath)
if d is None:
print("Error: Could not open image " + raster_filepath)
sys.exit(1)
global proj,geotrans,row,col
proj=d.GetProjection()
geotrans=d.GetGeoTransform()
row=d.RasterYSize
col=d.RasterXSize
inband=d.GetRasterBand(1)
Xcell_size=int(abs(geotrans[1]))
Ycell_size=int(abs(geotrans[5]))
cell_size = (Xcell_size+Ycell_size)/2
cell_dist=dist/cell_size
in_array[in_array == (inband.GetNoDataValue() or 0 or -999)]=0
out_array=np.zeros_like(in_array)
temp_array=np.zeros_like(in_array)
i,j,h,k=0,0,0,0
print("Running distance buffer...")
while(h<col):
k=0
while(k<row):
if(in_array[k][h]>=1):
i=h-cell_dist
while((i<cell_dist+h) and i<col):
j=k-cell_dist
while(j<(cell_dist+k) and j<row):
if(((i-h)**2+(j-k)**2)<=cell_dist**2):
if(temp_array[j][i]==0 or temp_array[j][i]>((i-h)**2+(j-k)**2)):
out_array[j][i]= in_array[k][h]
temp_array[j][i]=(i-h)**2+(j-k)**2
j+=1
i+=1
k+=1
h+=1
d,temp_array,in_array=None,None, None
return out_array

def export_array(in_array,output_path):
"""This function is used to produce output of array as a map."""
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(output_path,col,row,1)
outband=outdata.GetRasterBand(1)
outband.SetNoDataValue(np.nan)
outband.WriteArray(in_array)
# Georeference the image
outdata.SetGeoTransform(geotrans)
# Write projection information
outdata.SetProjection(proj)
outdata.FlushCache()
outdata = None

raster_buffer_array=raster_buffer("Input_//stream_order.tif",1000)
export_array(raster_buffer_array,"Output//buffer.tif")
print("Done")
• i get the following error when i try to use your code - File "../ran_alignRasters/code/python_rasterBuffer.py", line 36, in raster_buffer if(temp_array[j][i]==0 or temp_array[j][i]>((i-h)**2+(j-k)**2)): IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices . Any tip? Commented May 3, 2023 at 9:39