1

I need to solve problem very similar to Generating polygon representing rough 100km circle around latitude/longitude point using Python? but subtly different.

I need a shape that will be a circle once projected in Web Mercator - not in reality.

Working answer ( https://gis.stackexchange.com/a/268251/45061 ) is not generalizable, as it relies on https://github.com/jwass/geog that will not allow doing what I am trying to do.

3
  • 4
    A geometric object which is a circle when presented in Web Mercator is of very limited use ("How not to generate geometries" seminar is all I can think of). The easiest way to generate is a simple buffer in Web Mercator, then deproject. But first you need to select a GIS software stack.
    – Vince
    Mar 8, 2021 at 20:28
  • @Vince I know that it is fairly unusual. If anyone worries about using it for anything serious, it is for generating boundaries of laser-cut decoration plates on scales small enough that usual Web Mercator deficiencies are not relevant (city centers - so remembering to keep the same scale should be sufficient). Also, as it is for decorations most of Web Mercator issues are not really relevant. Mar 8, 2021 at 22:48
  • Web Mercator's deficiencies as a map projection are so numerous that I consider paper with ink jet ink an unconscionable waste of materials in map production. I can't imagine using any more exotic material being used. If the print cost is measured in dollars, a local Albers Equal Area or Lambert Conformal would be a more worthy projection.
    – Vince
    Mar 9, 2021 at 12:04

1 Answer 1

1
+200

Despite it being an abomination in so many ways, Pseudo/Spherical (Web) Mercator is a projection after all, so basic math on the Cartesian plane applies:

import math

R = 6378137.0


def LL2SM(ll):
  return [math.radians(ll[0]) * R, math.log(math.tan(math.pi / 4 + math.radians(ll[1]) / 2)) * R]

def SM2LL(xy):
  return [math.degrees(xy[0] / R), math.degrees(2 * math.atan(math.exp(xy[1] / R)) - math.pi / 2.0)]

def VERTEX(cxy, a, r):
  return [cxy[0] + r * math.cos(math.radians(a)), cxy[1] + r * math.sin(math.radians(a))]


def main():
  vertices = 32

  center = [13.0, 52.0]
  radius = 100000

  circle = []

  _angle = 360.0/vertices
  for vertex in range(0, vertices):
    sm = LL2SM(center)
    cv = VERTEX(sm, vertex*_angle, radius)
    ll = SM2LL(cv)

    circle.append(ll)

  # do sth. with your circle
  print(circle)


if __name__ == "__main__":
    main()

This POC script

  1. transforms a pair of Longitude/Latitude (center) into EPSG:3857 coordinates (LL2SM())
  2. generates vertices amount of points on a circle with a given radius (in meter) around center, using its parametric equations (VERTEX())
  3. transforms each vertex' coordinates back to Longitude/Latitude (SM2LL())

circle will then hold coordinate arrays ([Longitude, Latitude]) which, when transformed again to EPSG:3857, will form a circle around center.

Obviously, you can skip transforming back to Longitude/Latitude if you want.


You haven't specified any software environment other than Python, so I leave translating this into the framework of your choice to you; this includes generating an actual Polygon, for which there likely are built-in functions.

Edit:

...
# do sth. with your circle
wkt_polygon = ''.join(['POLYGON((', ','.join([' '.join(map(str, ll)) for ll in circle + [circle[0]]]), '))'])

print(wkt_polygon)

will output a Polygon as valid Well Known Text representation.

1
  • Note that this is obviously the basic and manual math version of the general workflow (project center -> create buffer -> deproject buffer), which should be covered in all spatial frameworks. Good to know the math, though. Keep it simple once in a while.
    – geozelot
    Mar 16, 2021 at 11:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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