# Can't get to plot a set of squares of defined length inside a circle plotted in EPSG 4326 using GeoPandas

I have been trying to plot a part of Quebec using WGS84 (EPSG:4326) degree coordinates.

I wanted to draw a circle within the map (e.g., 10 km) using the same coordinates and then subdivide the circle into squares of 1km side length. The resulting circle and squares look like this:

I got to plot the circle in the map using EPSG:4326; using the following code (As far as I understand, this implied converting it into an appropriate projection with meters as the coordinates and then converting it back to EPSG:4326, please correct me if my understanding is wrong as I am new to this):

``````df = pd.DataFrame(
{
'lat':[45.8],
'lon':[-74],
}
)
gdf = geopandas.GeoDataFrame(
df, geometry=geopandas.points_from_xy(df.lon, df.lat), crs="epsg:4326")
df_flat = gdf.to_crs('epsg:3857') #or 3797 or 3347
gdf['geom']=gdf_flat.geom.to_crs('epsg:4326')
``````

Now following this, I tried to do the same for the squares within the circle, taking the top right quadrant and the first square with one corner vertex at the center of the circle, I wanted to see if I could plot the squares using a similar method as described above.

I tried multiple methods, but they all led to different results, none of which helped me plot the square with one corner on the circle's origin.

Method 1 (using geopy):

``````import geopy as gp
d = distance.geodesic(kilometers=1, ellipsoid='WGS-84')
p1 = gp.Point((-74, 45.8))
p2 = d.destination(point=p1, bearing=0)
p3 = d.destination(point=p2, bearing=90)
p4 = d.destination(point=p3, bearing=180)
points = [(p.latitude, p.longitude) for p in [p1,p2,p3,p4]]
polygon_a = shapely.geometry.Polygon(points)
poly_df = pd.DataFrame()
poly_df['geom'] = pd.DataFrame.from_dict({"geom":[polygon_a]})
Polygdf_1 = geopandas.GeoDataFrame(poly_df, geometry='geom', crs='epsg:4326')
polygdf_flat = Polygdf_1.to_crs('epsg:6622') #6347
polygdf_flat['geom'] = polygdf_flat.geometry.buffer(1000, join_style=2)
Polygdf = polygdf_flat.to_crs('epsg:4326')
``````

Method 2 (Using centroid of square to create the square, I suspect I am doing this wrong, as I might need to set the CRS for p11)

``````import geopy as gp
p11 = gp.Point((-74, 45.8))
p21= distance.distance(kilometers=0.5).destination(p1, bearing=0)
p31= distance.distance(kilometers=0.5).destination(p2, bearing=90)
poly_df = pd.DataFrame()
poly_df['geom'] = pd.DataFrame.from_dict({"geom":[Point(p31[0], p31[1])]})
Polygdf_1 = geopandas.GeoDataFrame(poly_df, geometry='geom', crs='epsg:4326')
polygdf_flat = Polygdf_1.to_crs('epsg:3857') #6347 54032
polygdf_flat['geom'] = polygdf_flat.geometry.buffer(1000, cap_style=3)
Polygdf = polygdf_flat.to_crs('epsg:4326')**strong text**
``````

Method 3 (more like method 2b)

``````gdf_point = gdf.copy()
p111 = geopandas.GeoDataFrame(geometry=gdf_point['geometry'], crs="epsg:4326")
p111_mod = p111.to_crs('epsg:2035')#3347, 3348, 3573,32198
p211_mod = geopandas.GeoDataFrame(
p111_mod, geometry=geopandas.points_from_xy(p111_mod.geometry.x+500, p111_mod.geometry.y), crs="epsg:2035")
p311_mod = geopandas.GeoDataFrame(
p111_mod, geometry=geopandas.points_from_xy(p111_mod.geometry.x, p111_mod.geometry.y+500), crs="epsg:2035")
p411_mod = geopandas.GeoDataFrame(
p111_mod, geometry=geopandas.points_from_xy(p111_mod.geometry.x+500, p111_mod.geometry.y+500), crs="epsg:2035")
point_p411 = p411_mod.to_crs('epsg:4326')
p411gdf_m = geopandas.GeoDataFrame(columns=['geom'], geometry='geom')
p411gdf_m['geom'] = p411_mod.geometry.buffer(1000, cap_style=3)
p411gdf_llon = p411gdf_m.to_crs('epsg:4326')
``````

Upon plotting all of this on a map, the circle and the squares don't get plotted as desired in the first figure. (Green square is using Method 1, Blue Square is using Method 2, Thistle square is using Method 3, of the three, atleast the Blue square stays within the desired quadrant and does not move around so wildly).

I understand that the rotation is brought on by the differences in the type of projections that I have used, but I can't even seem to get them to line up to the right quadrant or consistently behave in the way I want. Can someone tell me what I am doing wrong?? I am sure I am missing some basic understanding of the mapping concept, but after many attempts, my brain is rather scrambled. If I can get this right for one or two squares (even if they look like rectangles on the map but stick to the right quadrant, I am sure I can figure it out for the rest of the squares.)

P.S: I did play around a bit with the different EPSG codes. This is the document I used for the EPSG projection codes: https://mern.gouv.qc.ca/wp-content/uploads/CO_codes_epsg_quebec.pdf

Let's start from the beggining. As a first step you should create a circle. I don't know why you created the "geom" column instead of using the "geometry" one provided by geopandas but in my case it turns out like this:

``````df = pd.DataFrame(
{
'lat':[45.8],
'lon':[-74],
}
)
gdf = gpd.GeoDataFrame(
df, geometry=gpd.points_from_xy(df.lon, df.lat), crs="epsg:4326")
gdf = gdf.to_crs('epsg:3857') #or 3797 or 3347
gdf=gdf.to_crs('epsg:4326')
gdf.plot()
``````

Now step 2 is to create a grid on top of this. I share here a piece of code that produce a grid based on an AOI (gdf). It will need to be adapted to solve projection issues:

``````from itertools import product

from shapely import geometry as sg
from shapely.ops import unary_union
import geopandas as gpd
import numpy as np
import ee

def set_grid(gdf, size):
"""
compute a grid around a given aoi (ee.FeatureCollection) that is fit for alert extraction
The grid cells are tailored to be adapted to always run without timeout

Args:
gdf: geodatframe of the aoi geometry
size: the square size in meters

Return:
the grid in 4326
"""

# project the gdf to 3857
proj_gdf = gdf.to_crs(3857)

# retreive the bounding box
aoi_bb = sg.box(*proj_gdf.total_bounds)
min_x, min_y, max_x, max_y = aoi_bb.bounds

# create numpy corrdinates table
lon = np.concatenate([np.arange(min_x, max_x, size), [max_x]])
lat = np.concatenate([np.arange(min_y, max_y, size), [max_y]])

# create the grid
squares = []
for ix, iy in product(range(len(lon) - 1), range(len(lat) - 1)):

# fill the grid values
square = sg.box(lon[ix], lat[iy], lon[ix + 1], lat[iy + 1])
squares.append(square)

# create a buffer grid in lat-long
grid = gpd.GeoDataFrame({"geometry": squares}, crs="EPSG:3857")
grid = grid.to_crs(4326)

return grid
``````

But this grid is covering the full with of the bounding box of the circle so we then need to filter it byt using the within method of geopandas:

``````grid = set_grid(gdf, 1000)

# filter outgeometry ouside of the bounds
geometry = gdf.dissolve().geometry.iloc[0]
grid = grid[grid.within(geometry)]
``````

finally the plotting of the result look like this on my machine:

``````from matplotlib import pyplot as plt

fig, ax = plt.subplots()
grid.plot(ax=ax, edgecolor='black')
gdf.plot(ax=ax, facecolor='none', edgecolor='red', linewidth=2)
``````

• Firstly, interesting approach, thanks for pointing me to the grid approach, that is the part that was helpful. I had not thought of using a shapely box as the bounds for EPSG:3857. That is a beautiful solution.. I am using "geom" and "geometry" in the 1st step only because this is an intermediate step in a larger optimization problem, which requires both the point and the polygon geometries and the same dataframe was being maintained for ease of use. I realize it is not the best way to manage the problem.. Jul 26, 2022 at 18:06
• The dissolve is not appropriate for our problem as we have a 2-D packing problem with slightly different constraints to decide on the number of squares themselves. That outputs 69 instead of 67 squares per quadrant for this example.. But I am marking this as the answer as it tackles the core of the issue with the grid.. Jul 26, 2022 at 18:07