3

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?

9

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

2
  • 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 '19 at 14:24
  • What can I do to buffer different distances for x and y direction? .buffer only takes one distance right? – till Kadabra Feb 21 '20 at 12:29

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.