Geographic Information Systems Asked by Hisham Sajid on February 18, 2021
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?
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.
First; fix your input lat/lon argument order.
Second; checkout your map projection (maps are not the truth, there is always a trade-off).
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:
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:
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
):
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:
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:
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...?:
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:
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.
Answered by s.k on February 18, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP