# Creating world grid with fix number of row and columns without overlap using shapely and geopandas

I'm trying to build a grid over the world and cut it on an AOI. The initial grid have non trivial properties.

It is a grid of squares in EPSG:4326. The dimension is not round, they simply cut the world horizontally into 2048 squares ( => 1024 vertically) or .17578125° per cell. The point [0,0] being the lower left corner of a cell.

I wrote a Python script that manages to create this kind of grid :

``````from shapely.geometry import Point
import geopandas as gpd
import decimal as d
from itertools import product
import numpy as np

# the size is based on the planet grid size
# the planet grid is composed of squared grid that split the world width in 2048 squares
diametre = 360/2048

# compute the longitudes and latitudes for the whole world
longitudes = np.linspace(-180, 180, 2048)
latitudes = np.linspace(-90, 90, 1024)

# filter with the geometry bounds
min_lon, min_lat, max_lon, max_lat = aoi_gdf.total_bounds
longitudes = longitudes[(longitudes > (min_lon - radius)) & (longitudes < max_lon + radius)]
latitudes = latitudes[(latitudes > (min_lat - radius)) & (latitudes < max_lat + radius)]

# create the grid
points = []
for i, coords in enumerate(product(longitudes, latitudes)):

x = d.Decimal(coords[0])
y = d.Decimal(coords[1])

points.append(Point(x, y))

# create a buffer grid in lat-long
grid = gpd.GeoDataFrame({'batch': batch, 'geometry':points}, crs='EPSG:4326') \
.buffer(d.Decimal(diametre)) \
.envelope \
.intersection(aoi_shp_proj)

# filter empty geometries
grid = grid[np.invert(grid.is_empty)]
``````

unfortunately when I watch my grid in QGIS I realized that the cells are actually overlaping. I also try without using the decimal lib but I end up with the same kind of problems (exact value of the Python float).

Is there a way to produce this kind of grid without creating ill shapes ?

## EDIT

To show what I mean by ill shape. Here is the grid produce on Singapore :

And here is a zoom on the corner

The tiles are in fact not aligned and overlap

• can you show us a picture of what the problem is? Commented Mar 3, 2021 at 17:01

Ok so I find a way to create the grid I wanted. 2 tricks where used:

• First do not create buffer but construct the squares with the points coordinates to make sure that the same value is used (we cannot do arythmetic with floats unfortunately)
• Second, in shp format the grid points where too precises and I was getting errors with the definitions of corners. GeoJson only use 13 decimals so the troncature was actually working. (reexport it to .shp afterward if needed)

NOTE: `aoi_shp` is a shapely geometry in EPSG:4326 projection. I added an extra step with a `aoi_gdf`to easily move from 1 projection to another

``````from pathlib import Path
from itertools import product

import geopandas as gpd
from shapely import geometry as sg
import numpy as np
from pyproj import CRS

# get the shape of the aoi in EPSG:4326 proj
aoi_gdf = gpd.GeoDataFrame({'geometry': [aoi_shp]}, crs="EPSG:4326")

# retreive the bb
aoi_bb = sg.box(*aoi_gdf.total_bounds)

# get the crs area of use
crs_4326 = CRS.from_epsg(4326)
crs_min_x, crs_min_y, crs_max_x, crs_max_y = crs_4326.area_of_use.bounds

# create the set of point
longitudes = np.linspace(crs_min_x, crs_max_x, 2048+1) # here you will have the points not the squares so +1
latitudes = np.linspace(crs_min_y, crs_max_y, 1024+1)

box_size = (crs_max_x-crs_min_x)/2048

# filter with the geometry bounds
min_lon, min_lat, max_lon, max_lat = aoi_gdf.total_bounds

# filter lon and lat
lon_filter = longitudes[(longitudes > (min_lon - box_size)) & (longitudes < max_lon + box_size)]
lat_filter = latitudes[(latitudes > (min_lat - box_size)) & (latitudes < max_lat + box_size)]

squares = []
for coords in product(range(len(lon_filter)-1), range(len(lat_filter)-1)):

ix = coords[0]
iy = coords[1]

squares.append(sg.box(
lon_filter[ix],
lat_filter[iy],
lon_filter[ix+1],
lat_filter[iy+1]
))

# create a buffer grid in 4326 from the drawn squares
grid = gpd.GeoDataFrame({'geometry':squares}, crs='EPSG:4326')

# cut the grid to the country extends