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?

Answer 1

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

Answer 2

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