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
Jul 31 '17 at 20:02
• Why the hell you guys down vote? his answer below is awesome! Jul 31 '17 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
Jul 31 '17 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. Aug 1 '17 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))
Ycell_size=int(abs(geotrans))
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")