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