8

I am using GDAL to convert a raster dataset into a numpy.array and then using numpy to buffer the data in the raster. For example, I need to know all the regions within 100ft of a school point feature and I need it represented as a raster.

So, I convert my points to a binary raster (school points are 1 and all other pixels are) using GDALRasterize in Python. Then, I open the raster dataset and use GDAL.ReadAsArray to create a numpy.array for the raster. Using scipy.ndimage.morphology.binary_dilation, I buffer all the non-zero pixels.

The issue is: some of my input datasets can be rather large and so I am running out of memory. The first thing that came to mind was to tile my input data and process individual tiles. I'm not sure how to handle any edge effects if I were to tile. What if a buffer spans the edge of a tile? The next tile wouldn't know to include those buffered cells.

Does anyone have any suggestions for how to handle these kinds of edge effects in numpy/GDAL?

3
  • 2
    Have you considered buffering the points and then rasterizing the buffer polygons? Example of buffering here.
    – user2856
    Oct 1, 2013 at 1:01
  • @Luke - I have considered buffering the vector data beforehand, but I don't think that is necessarily a viable option for me. I am trying to automate a process I have done previously manually. The input data can be points, lines, or polygons and can cover extents up to hundreds of sq miles. When I manually did this process I tried to buffer all building footprints in a 500 sq mile region. ArcGIS crashed trying to buffer so many polygons, but converting to raster worked just fine.
    – Brian
    Oct 1, 2013 at 14:35
  • 1
    Is the overall aim to create a distance-to-pixel layer? If so would using the gdal_proximity tool ('gdal.org/gdal_proximity.html') not be be better.
    – danclewley
    Nov 19, 2015 at 6:08

2 Answers 2

2

You could use XArray + Dask to split the procedure into chunks:

import numpy as np
import xarray as xr

# Create some dummy data.
array = np.zeros((11,11))
array[5,5] = 1
array[0,2] = 1
array[0,3] = 1

# Put the data in a xarray.DataArray and make it a dask chunked array.
da = xr.DataArray(array, dims = ["y", "x"]).chunk("auto")

# Create the buffers in two directions.
x_dir = da.rolling(x=3, center = True, min_periods = 2).sum()
y_dir = da.rolling(y=3, center = True, min_periods = 2).sum()

# Sum the buffers together.
test = (x_dir + y_dir) >= 1

# Make a plot to check the output.
img = test.plot()
img.axes.grid()

Buffered output

Instead of using gdal.ReadAsArray you'd have to open the data directly into XArray (doing something like da = xr.open_dataset(filepath, chunks = "auto")[var_name]), otherwise you are still loading the entire array into local memory.

0

You could consider making use of overlap in the tiles. Using an overlap of X pixels in each direction when reading in the data gives you more context, while writing out only the original tile (without overlap).

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.