6

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?

Improve this question

2 Answers 2

Reset to default
13

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

enter image description here

Improve this answer
4
  • 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! Commented 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? Commented 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. Commented Jun 21, 2022 at 9:29
  • @Irene, how did you plot the image, i have issue here. Thank you! Commented Dec 5, 2022 at 16:06
3

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

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.