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)
sq_buf = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 5, 5)

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?

2 Answers 2


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)
POLYGON ((13 15, 13 -5, -7 -5, -7 15, 13 15))

enter image description here

  • 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 5, 2019 at 14:24
  • What can I do to buffer different distances for x and y direction? .buffer only takes one distance right? Feb 21, 2020 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. Jun 21, 2022 at 9:29
  • @Irene, how did you plot the image, i have issue here. Thank you!
    – Lusia
    Dec 5, 2022 at 16:06

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

enter image description here

Resolution == 2

enter image description here

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.