2

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?

closed as off-topic by ahmadhanb, Jochen Schwarze, whyzar, Andre Silva, Vince Mar 5 at 19:15

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "When seeking help to debug/write/improve code always provide the desired behavior, a specific problem/error and the shortest code (as formatted text, not pictures) needed to reproduce it in the question body. Providing a clear problem statement and a code attempt helps others to help you." – ahmadhanb, Jochen Schwarze, whyzar, Andre Silva, Vince
If this question can be reworded to fit the rules in the help center, please edit the question.

4

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

  • 1
    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 at 14:24

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