# Define latitude longitude grid from center point

I would like to create a regular squared grid from a particular point. I have a polygon, I take the center of this polygon and I would like to create a grid that contains this polygon using Python.

I make the problem simpler by just working on a single point that is representing the center of my polygon.

``````import shapely
import numpy as np
import geopandas as gpd

lat, lon = 50, -10 # my point
``````

I'm using the Haversine distance to calculate the distance between the 4 points of each mesh.

``````from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles.
return c * r
``````

Following this case: How can I create grid of coordinates from a center point in Python?

I'm doing the same but i'm getting rectangular grid and not a squared one.

``````dist = 7000 # x and y size of each mesh in meters
coors = 2 #number of mesh in each direction

#Creating the offset grid
mini, maxi = -dist*coors, dist*coors
n_coord = coors*2+1
axis = np.linspace(mini, maxi, n_coord)
X, Y = np.meshgrid(axis, axis)

#avation formulate for offsetting the latlong by offset matrices
dLat = X/R
dLon = Y/(R*np.cos(np.pi*lat/180))
latO = lat + dLat * 180/np.pi
lonO = lon + dLon * 180/np.pi

#stack x and y latlongs and get (lat,long) format
output = np.stack([latO, lonO]).transpose(1,2,0)
list_polygon = []
len_x, len_y, _ = output.shape

for x_index in range(len_x - 1):
for y_index in range(len_y - 1):
p1 = output[x_index, y_index, :]
p2 = output[x_index, y_index+1, :]
p3 = output[x_index+1, y_index+1, :]
p4 = output[x_index+1, y_index, :]
print(haversine(*p1, *p2))
print(haversine(*p2, *p3)) # this distance is not equal to 7000 meters
print(haversine(*p3, *p4)) # and this distance too
print(haversine(*p4, *p1))
print()
polygon = shapely.geometry.Polygon([p1, p2, p3, p4])
list_polygon.append(polygon)
gdf_grid = gpd.GeoDataFrame({'geometry' : list_polygon}, index =
range(len(list_polygon)), crs = 'EPSG:4326')
``````

the Output is a rectangular grid, as expected with the distance.

``````import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1, figsize = (6, 6))
gdf_grid.plot(facecolor = 'none', ax = ax)
ax.axis('equal')
``````

I guess the issue comes from dLon but I cant figure why. Do you have any idea?

Reproject to a projected coordinate system. Then create a point grid, square buffer and go back to lat/long:

``````import geopandas as gpd
import shapely
import numpy as np

lat, lon = 50, -10 # my point
pnt = shapely.geometry.Point(lon, lat)

df = gpd.GeoDataFrame(geometry=[pnt], crs=4326)
df = df.to_crs(df.estimate_utm_crs())
dist = 7000
coors = 2

xcenter = df.geometry.iloc[0].x
ycenter = df.geometry.iloc[0].y

xcoords = np.arange(start=xcenter-dist*(coors-1), stop=xcenter+dist*coors, step=dist)
ycoords = np.arange(start=ycenter-dist*(coors-1), stop=ycenter+dist*coors, step=dist)

coords = np.array(np.meshgrid(xcoords, ycoords)).T.reshape(-1,2)

centerpoints = gpd.points_from_xy(x=coords[:,0], y=coords[:,1])
squares = [p.buffer(distance=dist, cap_style=3) for p in centerpoints]
df2 = gpd.GeoDataFrame(geometry=squares, crs=df.crs)

df = df.to_crs(4326)
df2 = df2.to_crs(4326)
ax = df.plot(color="red", markersize=200, figsize=(15, 15), zorder=1)
df2.boundary.plot(ax=ax, zorder=0)
``````