-1

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?

  • 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 Jul 31 '17 at 20:02
  • Why the hell you guys down vote? his answer below is awesome! – nickves Jul 31 '17 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 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. – Aakash Bikram Rana Aug 1 '17 at 4:17
4

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")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.