How to create a squared buffer around a point with Python GDAL?

Question

I have a large collection of points with X and Y coordinates in a projected coordinate reference system (EPSG:28992). I would like to compute a square buffer around it (say, 500 meters around the point), so I can crop this polygon with a raster layer (in TIFF format) and proceed with a further analysis with the values within the cropped windows.

I am aware that this is possible in Python for QGIS, but for the sake of the workflow I am implementing, I prefer to have it entirely in the Python GDAL environment, so it becomes automatic. Since I know the X/Y coordinates and the size of the grid cells in the raster, I could add/subtract to the coordinates the desired buffer, but I wonder whether there is a less handcrafted version of this "square-buffer" operation, which is less prone to errors. Other Python libraries (e.g. shapely) would be fine for me as well.

Note that a band.GetRasterBand(1).ReadAsArray(X, Y, 5, 5) does not return the window around the point, but the window leaving the point in a corner.

For instance, the following code:

path_curr_lyr = r"path_to_tif_file"
tif = gdal.Open(path_curr_lyr)
dat = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 1, 1)
print(dat)
sq_buf = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 5, 5)
print(sq_buf)

Accesses the raster grid cell X/Y (9550, 9620) and returning the value of forest coverage in that pixel, and the value of a window of width= 5 pixels.

>>> [[44]]
>>> [[ 44  87  99  86  43]
     [ 96  99  99  79  14]
     [ 99  99  98  34   9]
     [ 99  86  27  93  71]
     [ 98  98 100  97  65]]

As you can see, this approach returns a window in which the "44" is not a central value, but one in the corner. I would like to have a procedure in which a window around the given grid cell coordinate is created.

Any advice on how to implement this?

score 12 · Accepted Answer · answered Mar 05 '19 at 12:36

12

With shapely, you can use the styles of caps of buffer

from shapely.geometry import Point
test = Point(3,5)
# buffer with CAP_STYLE = 3
buf = test.buffer(10, cap_style=3)
print(buf.wkt)
POLYGON ((13 15, 13 -5, -7 -5, -7 15, 13 15))

answered Mar 05 '19 at 12:36

gene

54,868
3
110
187

2

This is is a really really nice solution. After the buffer calculation, I opened the TIF file with rasterio and then I used rasterio.mask.mask() function to crop the raster layer with the freshly calculated squared window. It worked perfectly. Concretely, I added out_img, out_transf = mask(tif, [buf], crop=True) in case someone is interested. Thanks! – Irene Mar 05 '19 at 14:24
What can I do to buffer different distances for x and y direction? .buffer only takes one distance right? – till Kadabra Feb 21 '20 at 12:29
1

@tillKadabra there is a single_sided parameter you can use for linestrings/polygons (use postitive/negative values for buffer to indicate side to buffer), but I this will not work for Points. – ryanjdillon Jun 21 '22 at 09:29
@Irene, how did you plot the image, i have issue here. Thank you! – Lusia Dec 05 '22 at 16:06

score 1 · Answer 2 · answered Jun 21 '22 at 09:37

You can also set the resolution equal to 1 to result in a squared buffer around a Point.

Interestingly, you get a rotated square with this method, so this could be helpful if that is what you want.

Borrowing the example code from this answer, but using the resolution parameter instead, I get the following:

import geopandas as gpd
from shapely.geometry import Point
import matplotlib.pyplot as plt
Generate some sample data
p1 = Point((1,2))
p2 = Point((6,8))
points = gpd.GeoSeries([p1,p2])
Buffer the points using a resolution of 1
buffer = points.buffer(2, resolution=1)
Plot the results
fig, ax1 = plt.subplots()
buffer.boundary.plot(ax=ax1, color = 'slategrey')
points.plot(ax = ax1, color = 'red')

Resolution == 1

Resolution == 2

How to create a squared buffer around a point with Python GDAL?

2 Answers2

Generate some sample data

Buffer the points using a resolution of 1

Plot the results

Linked

Related