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

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

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