3

I'm trying to create a circle around a point using the code from this answer here. Adding below for better visibility as well

from functools import partial
import pyproj
from shapely.ops import transform
from shapely.geometry import Point

proj_wgs84 = pyproj.Proj('+proj=longlat +datum=WGS84')


def geodesic_point_buffer(lat, lon, km):
    # Azimuthal equidistant projection
    aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
    project = partial(
        pyproj.transform,
        pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
        proj_wgs84)
    buf = Point(0, 0).buffer(km * 1000)  # distance in metres
    return transform(project, buf).exterior.coords[:]

# Example
b = geodesic_point_buffer(45.4, -75.7, 100.0)

print(b)
# [(-74.42290765358695, 45.39286001598599),
#  (-74.43102886629593, 45.304749544147974),
#  ...
# (-74.42290765358695, 45.392860015985995),
# (-74.42290765358695, 45.39286001598599)]

However when I try plot it on a map it looks more like an Ellipse and less like a circle. The coordinates I am trying to create the buffer around are (31.1482371,72.6859082). Is this some EPSG issue?

enter image description here

8
  • 3
    31.1482371E 72.6859082N is north of Finland in the Barents Sea and a buffer about it would have a significant egg-shape; 72.6859082E 31.1482371N is in Pakistan, and would have a less-pronounced wobble. Neither would be a proper ellipse. – Vince Nov 26 '20 at 23:03
  • 2
    I suspect your latitude and longitude coordinates are passed in the wrong order somewhere. Can you post the exact code that you used to generate this, with your actual coordinates? – FSimardGIS Nov 27 '20 at 3:08
  • Note that 72.6859082 is your longitude and 31.1482371 is your latitude. Make sure you pass those values in the correct order in your function. The function call should be in your case : geodesic_point_buffer(31.1482371, 72.6859082, 100.0) so latitude first, then longitude, and the buffer radius you need. – FSimardGIS Nov 27 '20 at 5:51
  • @FSimardGIS that's exactly how I'm running my code. I still suspect it's a projection issue. Probably ` aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0' ` something I need to change here but I can't figure it out. – Hisham Sajid Nov 27 '20 at 7:25
  • 1
    A very nice online tool to visualize Tissot's indicatrix as well as geodesic buffers/distance circles (on a global level) as well as different projections: schweizerweltatlas.ch/swatools/MapProjector/… – Babel Jan 19 at 21:55
5

The Earth is actually not flat.

It is not a projection issue, things are how they are.

Let's try to understand what's happening.

TL;DR:

First; fix your input lat/lon argument order.
Second; checkout your map projection (maps are not the truth, there is always a trade-off).

Details

Here is a revised piece of code:

import numpy as np
from functools import partial
import matplotlib.pyplot as plt
import pyproj
from pyproj.crs import ProjectedCRS
from pyproj.crs.coordinate_operation import AzumuthalEquidistantConversion
from shapely.ops import transform
from shapely.geometry import Point
import geopandas as gpd

def geodesic_point_buffer(lon, lat, km):
    proj_crs = ProjectedCRS(
        conversion = AzumuthalEquidistantConversion(lat, lon)
    )
    proj_wgs84 = pyproj.Proj('EPSG:4326')
    Trans = pyproj.Transformer.from_proj(
        proj_crs,
        proj_wgs84,
        always_xy=True
    ).transform

    return transform(Trans, Point(0, 0).buffer(km * 1000))

Then if you call your function for lon=0, lat=70 and km=100:

geobuf = geodesic_point_buffer(0, 70, 100)
gpd.GeoSeries(geobuf).plot()

it will show something like:

Result

You can even try to figure out what's happening on a representative portion of the Earth:

km = 100 
longitudes = np.arange(0, 90, 30)
latitudes = np.arange(0, 90, 30)

cases = np.flip([(x,y) for x in longitudes for y in latitudes])

figsize = (10, 10)
fig, axs = plt.subplots(
    len(longitudes),
    len(latitudes),
    figsize=figsize
)
axs = axs.flatten()

for ax, case in zip(axs, cases):
    lon, lat = case
    geobuf = geodesic_point_buffer(float(lon), float(lat), km)
    g = gpd.GeoSeries(geobuf)
    g.plot(ax=ax, aspect='equal', color=np.random.rand(3,), alpha=0.2)
    ax.tick_params(labelbottom=False, labelleft=False)
    ax.set_title('lon:{}, lat:{}'.format(lon, lat))

And here's the result:

Results of the loop

Argh, the northern we go, the more ellipsoid shapes we have. What's wrong ?

Let's actually map one of those 100 km radius "circle" (let's hope these ellipsoidal shapes are just an illusion) into QGIS.
Let's pick the one at lon=0 lat=80 (km=100):

The culprit

Here the WKT string if you want to compare:

POLYGON ((5.1429335111222647 79.9603916313436685, 5.0744480577032895 79.8733291528321701, 4.9593109530228228 79.7878284927520127, 4.7996980921410275 79.7046644066590346, 4.5980578200420545 79.6245727688073117, 4.3570640056511660 79.5482464155426641, 4.0795745337711660 79.4763316883885125, 3.7685952499268129 79.4094255861472789, 3.4272491771699012 79.3480734400190357, 3.0587506714144381 79.2927670327295289, 2.6663840815509183 79.2439430909645495, 2.2534864210693830 79.2019820892817989, 1.8234335286860917 79.1672073125885447, 1.3796291873385931 79.1398841328943945, 0.9254966762295718 79.1202194641839469, 0.4644722432153828 79.1083613668266992, 0.0000000000000070 79.1043987799610875, -0.4644722432153690 79.1083613668266850, -0.9254966762295578 79.1202194641839469, -1.3796291873385780 79.1398841328943945, -1.8234335286860803 79.1672073125885447, -2.2534864210693617 79.2019820892817989, -2.6663840815509010 79.2439430909645353, -3.0587506714144221 79.2927670327295289, -3.4272491771698941 79.3480734400190357, -3.7685952499267992 79.4094255861472789, -4.0795745337711589 79.4763316883885125, -4.3570640056511607 79.5482464155426641, -4.5980578200420510 79.6245727688073117, -4.7996980921410204 79.7046644066590346, -4.9593109530228228 79.7878284927519985, -5.0744480577032878 79.8733291528321701, -5.1429335111222647 79.9603916313436685, -5.1629158879165384 80.0482072395145821, -5.1329246533944488 80.1359391827525513, -5.0519298456202693 80.2227293444878882, -4.9194033618640605 80.3077060844526187, -4.7353796297580724 80.3899930802279954, -4.5005128730321271 80.4687192006217202, -4.2161276595591222 80.5430293477546684, -3.8842590173902676 80.6120961425510671, -3.5076782050576112 80.6751322581024226, -3.0899003094051007 80.7314031313892855, -2.6351702884681485 80.7802397121709390, -2.1484249195528449 80.8210508460416719, -1.6352293481534619 80.8533348449872733, -1.1016884971873431 80.8766897811666183, -0.5543353614832844 80.8908220542790701, -0.0000000000000161 80.8955528328694697, 0.5543353614832505 80.8908220542790701, 1.1016884971873147 80.8766897811666183, 1.6352293481534266 80.8533348449872733, 2.1484249195528125 80.8210508460416719, 2.6351702884681085 80.7802397121709390, 3.0899003094050776 80.7314031313892997, 3.5076782050575912 80.6751322581024226, 3.8842590173902507 80.6120961425510671, 4.2161276595591026 80.5430293477546684, 4.5005128730321156 80.4687192006217202, 4.7353796297580564 80.3899930802280096, 4.9194033618640498 80.3077060844526187, 5.0519298456202595 80.2227293444879024, 5.1329246533944417 80.1359391827525798, 5.1629158879165367 80.0482072395145821, 5.1429335111222700 79.9603916313436827, 5.1429335111222647 79.9603916313436685))

And here is how it looks like when loaded into QGIS:

Shape imported into QGIS

Okay... it's definitely an ellipse... or? :/

Wait, watch this amazing video: https://www.youtube.com/watch?v=YES0UtBelNM

And remember, WGS84 is a Geographic 2D projection according to the EPSG authority itself:

EPSG:4326

Source: https://epsg.org/crs_4326/WGS-84.html?sessionkey=m6fd9wl8om

What does all this mean? It means that you actually correctly measured a true circle in a 3D dimensional space, onto a "sphere" - the Earth (which actually isn't really a sphere, you know...) -, and that you tried to represent it on a flat surface (your map), which has introduced some de-facto distortions because we cannot represent a sphere onto a flat plane. Hence you need to take a projected coordinate system to check your shape, such as the well know Web Mercator, also known as "WGS84 / Pseudo-Mercator" or EPSG:3857. And guess what happens if we switch to this projection into QGIS...?:

Yeah circle !

Bingo, it became a circle!
An other way to convince yourself that it's a circle even when you have the WGS84 Geographic projection active on your QGIS project is to measure distances:

Circle in WGS84

Here you can see the circle on the first previous WGS84 map. Notice the distance measurement in kilometers in the box when moving the cursor. They are measured in ellipsoidal default metric as explained here: https://docs.qgis.org/3.16/en/docs/user_manual/introduction/general_tools.html#measuring
Simply imagine you are measuring true, geodetic, distances (as you wanted) on the surface of the true Earth, where your circle actually lies. It's what it means.

Hope this helps to make it clear.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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