I'm new to GDAL Python, and I'm working on a little recursive algorithm which requires finding the neighbors of a a pixel on a DEM raster. I've been able to do this by creating a little PixelLite class which is essentially just a tuple representing the raster coordinates of a pixel. Here's a bit of its implementation (feel free to skip):

    # a lighweight version of a pixel for quick computation 
    class PixelLite:

        def __init__(self, x, y, raster_meta):
            # set x and y
            self.x = x
            self.y = y

            # set max_x and max_y
            self.raster_meta = raster_meta
            self.max_x = raster_meta.x_max
            self.max_y = raster_meta.y_max

            # set error flag
            self.flag = False

            # the set of all neighbors
            self.neighbors = set()

            if self.x > self.max_x | self.x < 0 | self.y > self.max_y | self.y < 0: # something has gone terribly wrong
                self.flag = True

            # !!! Below, I'm going to set out of bounds values to the value of the pixel itself, which will later be removed from the set !!!

            # neighbor columns
            if self.x > 0: # if left value will not be negative
                left = self.x - 1 # add legit left
                left = self.x # add left value that will be removed

            center = self.x

            if self.x < self.max_x: # if right value will not be larger than raster dimensions
                right = self.x + 1 # add legit right value
                right = self.x # add right value that will be removed

            # Note: bottom and top are arbitrary names, what matters is that values do no exceed bounds; it may be the case that, because of raster indexing 'bottom' actually represents the pixel above the current pixel
            # neighbor columns
            if self.y > 0: # if bottom value will not be negative
                bottom = self.y - 1 # add legit bottom value
                bottom = self.y # add bottom value that will be removed

            middle = self.y

            if self.y < self.max_y: # if top value will not exceed raster extent
                top = self.y + 1 # add legit top value
                top = self.y # add top value that will be removed

            # top row of neighbors
            self.neighbors.add((left, top))
            self.neighbors.add((center, top))
            self.neighbors.add((right, top))

            # middle row of neighbors
            self.neighbors.add((left, middle))
            self.neighbors.add((right, middle))

            # bottom row of neighbors
            self.neighbors.add((left, bottom))
            self.neighbors.add((center, bottom))
            self.neighbors.add((right, bottom))

            # the method I'm using will sometimes add the center, middle element as I've used that as a default to set values to if there is an out of bounds error. So let's just pop that off...; a pixel cannot be its own neighbor
            self.neighbors.discard((center, middle))


Basically, what this does is create the tuples representing each neighbor for a given pixel whenever the PixelLite class is instantiated. Although this is pretty lightweight, because it is implemented in very basic Python without any kind of binding, I would imagine it lacks the optimization of a more sophisticated solution such as I would expect to find in GDAL.

Is there a faster getPixelNeighbors function provided by GDAL or another raster library i.e. one which would allow me to do the same but which might have some neat little optimization tricks up its sleeve?

  • There may be optimised neighbourhood code buried in the GDAL C++ lbrary for use by other function, but I'm not aware of it. I would use the gdal Dataset or RasterBand ReadAsArray to read as a numpy array and use the power of numpy indexing. Even easier, use rasterio, your open dataset object has a read method for reading as a numpy array and an index method to transform map coordinates to pixel indexes which you can use with your numpy array. – user2856 Apr 20 at 3:05
  • @user2856 Thanks! I'm currently using a gdal_array which is purportedly an easier way of interfacing between numpy and gdal. But, if I understand you correctly, it sounds like I might still gain some efficiency by dealing with a numpy array directly. Is that correct? – LIAM M Apr 20 at 3:21

Your Answer

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

Browse other questions tagged or ask your own question.