Creating square polygon around center of the point using Python

Question

How do I create a square polygon using center of the point (lat/lon) 25 kilometers (KM) by sides?

I am not looking for a radius i.e. buffer. Here I am creating edges of 25 KM from the center point my_lat and my_lon. I may be coming close to it, yet I don't know enough GeoPandas functions that can calculate the square/rectangle polygons. Any suggestions?

I have provided a picture below for a better understanding of the questions. Black lines are currently 35 KM apart.

import geopy
import geopy.distance as distance
import plotly.graph_objects as go
from haversine import haversine, Unit
from shapely.geometry import Polygon
Each side is 25KM appart
d = distance.distance(kilometers=25) # 25 KM
print(d)
Going clockwise, from lower-left to upper-left, upper-right...
my_lat = 38.696211
my_lon = -120.5
center_point = geopy.Point((my_lat,my_lon))
p2 = d.destination(point=center_point, bearing=45)
p3 = d.destination(point=center_point, bearing=135)
p4 = d.destination(point=center_point, bearing=-135)
p5 = d.destination(point=center_point, bearing=-45)
print('p2','-->',p2)
print('p3','-->',p3)
print('p4','-->',p4)
print('p5','-->',p5)
points = [(p.latitude, p.longitude) for p in [p2,p3,p4,p5]]
polygon = Polygon(points)[![enter image description here][1]][1]

Answer 1

Maybe you can adapt from this code which makes use of Shapely's buffer method through GeoPandas with a squared caps style:

import geopandas as gpd
from math import sqrt
from shapely import wkt
def square_poly(lat, lon, distance=25000/sqrt(2)):
    gs = gpd.GeoSeries(wkt.loads(f'POINT ({lon} {lat})'))
    gdf = gpd.GeoDataFrame(geometry=gs)
    gdf.crs='EPSG:4326'
    gdf = gdf.to_crs('EPSG:3857')
    res = gdf.buffer(
        distance=distance,
        cap_style=3,
    )
return res.to_crs('EPSG:4326').iloc[0]


my_lat = 38.696211
my_lon = -120.5
geom = square_poly(my_lat, my_lon, distance=25000/sqrt(2))
print(geom.wkt)

Results in the following ~35'355.34m sided square (50km diagonal):

POLYGON ((
    -120.34119879273891 38.820043441823955,
    -120.34119879273891 38.572163797794275,
    -120.65880120726112 38.572163797794275,
    -120.65880120726112 38.820043441823955,
    -120.34119879273891 38.820043441823955
))

Shown in QGIS:

EDIT:

The diagonal of a square with a side of 1 is sqrt(2) by definition. Hence, if the diagonal is 2·25'000, the side is therefore (50'000/2)/sqrt(2) or equivalently 25'000/sqrt(2).

From your initial drawing, it's not possible to draw a square with a side of 25 and a diagonal of 50 (or 25, it's not 100% clear from your drawing because the dashed red line is not the actual diagonal of the square but half of it).

So I let you figure out what you want to achieve by replacing one or the other element in the rule of three. Also, you may want to transform your input WGS84 coordinates to some local, more suited, Cartesian reference system. The EPSG:3857 was just used as a "default" placeholder here because it's the default CRS for many web maps and it may be enough depending on your goal.