0

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?

4
  • 1
    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!
    – nickves
    Commented Jul 31, 2017 at 20:36
  • 2
    @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
  • 1
    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

1 Answer 1

6

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)
     in_array = inband.ReadAsArray(0,0,col,row).astype(int)
     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")
1
  • 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?
    – Momchill
    Commented May 3, 2023 at 9:39

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.