GeoPandas – Plotting Squares of Defined Length Inside a Circle Plotted in EPSG 4326

geopandaspythonwgs84

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],
        'rad':[10000]
    }
)
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)
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

Best Answer

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],
        'rad':[10000]
    }
)
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)
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)

Related Solutions

GeoPandas – Fixing Warnings by Correctly Using EPSG Codes

You get that warning because EPSG:4326 (WGS84) is geographic CRS and its unit is degree. The key message here is "Geometry is in a geographic CRS". You don't get a warning when using EPSG:4328 because it is not geographic, but it is geocentric and its unit is meter.

When you use area, length or distance, GeoPandas checks CRS before calculation and warns if the planar operation is done in a geographic CRS.

Why checking?
One of the packages that Geopandas depends on is Shapely, where geometric calculations are made. Shapely uses 2D Cartesian coordinates (x, y) and calculations are performed in 2D plane. So Shapely doesn't care if the coordinates are in degree or meter. They are just numbers to Shapely. Therefore, in the table, the unit of calculated_area1 is square degree (no meaning), and the other is square meter.

GeoPandas Attribute Sum – Get the Sum of An Attribute Within a Set Radius for Each Row in GeoPandas

You can buffer your cells to create the coverage areas, then use a spatial join to link each coverage area with the cells it contains.

import geopandas
from shapely.geometry import Point

# create example data
cells = geopandas.GeoDataFrame({
    'index': [0,1,2,3],
    'pop': [1,2,3,4],
    'geometry': [
        Point([0,0]),
        Point([0,1]),
        Point([1,0]),
        Point([1,1])
    ]
}).set_index('index')

# create "coverage" buffers
coverage = cells[['geometry']].copy()
# take care here - your data should be in a projected CRS
# if you want to interpret this as a distance in meters
coverage['geometry'] = coverage.buffer(1.1) 

# join cells with coverage, then groupby and sum
cell_coverage_pop = cells.sjoin(coverage, predicate='within') \
    .groupby('index_right') \
    .sum()
# rename index and column name
cell_coverage_pop.index.rename('index', inplace=True)
cell_coverage_pop.rename(columns={'pop': 'coverage'}, inplace=True)
# cell_coverage_pop has "index" index and "coverage" 
# with population total

# join back to original cells if desired
cells = cells.join(cell_coverage_pop)
# cells now has "index" index, and original "pop", original 
# "geometry" and total "coverage"

Best Answer

Related Solutions

GeoPandas – Fixing Warnings by Correctly Using EPSG Codes

GeoPandas Attribute Sum – Get the Sum of An Attribute Within a Set Radius for Each Row in GeoPandas

Related Question