2

Is there an equivalent to the very good st_make_grid method of the sf package from r-spatial in python? The method create rectangular grid geometry over the bounding box of a polygon.

I would like to do exactly the same as the solution proposed in this question, e.g. divide a polygon into several squares of the same area that I choose. Thanks for your help.

Alternatively, I could use rpy2 to run a script in r that executes the st_make_grid method which takes a shapely polygon as input and outputs the square polygons, to be read with shapely. Would this be effective on many polygons to process?

Tim
  • 513
  • 5
  • 20
  • 1
    This can be as efficient in Python when using rpy2 as in R if you avoid frequently moving objects bidirectionally between R and Python. – krassowski Aug 13 '21 at 14:13
  • If have large arrays or data table Apache Arrow may provide very substantial performance gains when sharing data between Python and R. Check the rpy2 extension rpy2-arrow. – lgautier Aug 14 '21 at 20:38

1 Answers1

3

Would this be effective on many polygons to process?

Certainly not. There's no built-in Python version but the function below does the trick. If you need performance, make sure that you have pygeos installed in your environment.

def make_grid(polygon, edge_size):
    """
    polygon : shapely.geometry
    edge_size : length of the grid cell
    """
    from itertools import product
    import numpy as np
    import geopandas as gpd
    
    bounds = polygon.bounds
    x_coords = np.arange(bounds[0] + edge_size/2, bounds[2], edge_size)
    y_coords = np.arange(bounds[1] + edge_size/2, bounds[3], edge_size)
    combinations = np.array(list(product(x_coords, y_coords)))
    squares = gpd.points_from_xy(combinations[:, 0], combinations[:, 1]).buffer(edge_size / 2, cap_style=3)
    return gpd.GeoSeries(squares[squares.intersects(polygon)])
martinfleis
  • 7,124
  • 2
  • 22
  • 30
  • Firstly, thank you very much for your solution, it is very succinct and has helped me learn new things. Secondly, the area is the square root of the side length, right? For example with the following polygon: `POLYGON ((39.31869506835937 25.14777088272356, 39.88861083984375 25.14777088272356, 39.88861083984375 25.50526349022746, 39.31869506835937 25.50526349022746, 39.31869506835937 25.14777088272356))` an `edge_size` length greater than 0.8 returns an error. A length of 0.01 gives me an area of 1.36 km2, and 0.02 an area of 5.45km2. **How to choose the length for an area**? As 8km2 ? – Tim Aug 13 '21 at 23:19
  • For the calculation of the area I proceeded as follows: `result = make_grid(polygon, 0.02) \ gdf = gpd.GeoDataFrame(geometry=result) \ gdf.crs = "epsg:4326" \ gdf = gdf.to_crs(epsg=3395) \ print(gdf.area*1e-6) # in km2` – Tim Aug 13 '21 at 23:27
  • The function above assumes projected CRS, it will not work properly if your geometries are in EPSG:4326. – martinfleis Aug 14 '21 at 11:12
  • Indeed the calculations are in meters.. thanks again very much ! – Tim Aug 14 '21 at 11:31
  • Hi @martinfleis, I am facing a problem. I have an input polygon in 4326 projection. I project it in 3395 and execute the `make_grid` function. I then re-project the square polygons obtained in 4326. (All in one projection as you explained [here](https://stackoverflow.com/a/69031632/14864907)). Then I calculate the center points of each square with `.centroid` but even here they are inaccurate. Is it possible ? – Tim Sep 05 '21 at 10:46
  • A case where it works: from a point in 4326 if we create the square of this point in 3395 and re-project this square in 4326. The center of the square is equal to the starting point. But we created the square with its center in 4326 at the beginning. In the case of `make_grid` the squares are generated without initial center points in 4326. – Tim Sep 05 '21 at 11:01