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:

enter image description here

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(
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_flat['geom'] = gdf_flat.geometry.buffer(df.rad)

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). enter image description here

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

1 Answer 1


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(
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['geometry'] = gdf.geometry.buffer(df.rad)

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
        gdf: geodatframe of the aoi geometry
        size: the square size in meters
        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])

    # 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)

enter image description here

  • 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.