0

I'm trying to draw circles around lat-long (WGS84) points by:

  • transforming the point into AEQD space
  • creating a shapely buffer on the projected point
  • projecting the shapely buffer back to WGS space

```

import pyproj
from shapely.ops import transform as shapely_transform
from shapely.geometry import Point

from functools import partial

def non_euclidean_circle(center, radius):
    lat, lon = center
    # proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj='aeqd', lat_0=lat, lon_0=lon, x_0=0, y_0=0)
    WGS84 = pyproj.Proj(init='epsg:4326')

    # transform the given lat-long onto the flat AEQD plane
    tx_lat, tx_lon = pyproj.transform(WGS84, AEQD, lat, lon)
    circle = Point(tx_lat, tx_lon).buffer(radius)

    def inverse_tx(x, y, z=None):
      return pyproj.transform(AEQD, WGS84, x, y)

    # inverse projection from AEQD to EPSG4326-WGS84
    return shapely_transform(inverse_tx, circle)

c = non_euclidean_circle((42, -72), 100)
print(c, c.area)

When I plot the polygon c on a map, it looks like: elongated circle

Complete plotting code: https://codepen.io/anon/pen/yEoJrG?editors=0010

My question is, is this correct? Are these points on the manifold of the Earth actually equidistant from the center, or am I doing something wrong here?

closed as off-topic by Vince, BERA, Andre Silva, nmtoken, neogeomat Jun 17 '18 at 13:32

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This problem cannot or can no longer be reproduced. Changes to the system or to the asker's circumstances may have rendered the question obsolete, or the question does not include a procedure to enable potential answerers to reproduce the same symptoms. Such questions are off-topic as they are unlikely to help future readers, but editing them to include more details can lead to re-opening." – Vince, BERA, Andre Silva, nmtoken, neogeomat
If this question can be reworded to fit the rules in the help center, please edit the question.

  • Please Edit the question to specify the exact input paramaters. It's possible a far-north spheroidal circle would look like that, but I think you did something wrong. – Vince Jun 14 '18 at 13:37
  • Note: if you leave x_0,y_0 equal to zeroes, the center of your AEQD is always 0,0. So just build the buffer then unproject it using the customized AEQD definition.You don't need to project the original lat/lon values to AEQD. Also what ellipsoid/sphere is AEQD using? – mkennedy Jun 15 '18 at 0:38
  • @Vince I have specified the exact params (42, -72 or 42°00'00.0"N 72°00'00.0"W is the lat-long somewhere near MA, USA). – ixaxaar Jun 15 '18 at 7:18
  • 1
    No, {42,-72} is near Antarctica (latitude is a Y value) – Vince Jun 15 '18 at 9:53
  • @Vince thank you so much, I think THAT is exactly what confused me! Would you like to post it as an answer? – ixaxaar Jun 16 '18 at 11:06
0

Had the order of lat-long wrong. :/

This works as expected:

import pyproj
from shapely.ops import transform as shapely_transform
from shapely.geometry import Point

from functools import partial


def non_euclidean_circle(center, radius):
    lat, lon = center
    # proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj='aeqd', lat_0=lat, lon_0=lon, x_0=lon, y_0=lat)
    WGS84 = pyproj.Proj(init='epsg:4326')

    # transform the given lat-long onto the flat AEQD plane
    tx_lon, tx_lat = pyproj.transform(WGS84, AEQD, lon, lat)
    circle = Point(tx_lat, tx_lon).buffer(radius)

    def inverse_tx(x, y, z=None):
      y, x = pyproj.transform(AEQD, WGS84, y, x)
      return (x, y)

    # inverse projection from AEQD to EPSG4326-WGS84
    return shapely_transform(inverse_tx, circle)

c = non_euclidean_circle((42, -72), 100)

as expected

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