Bug for lat near 90 or -90 and lon near 180 or -180

[anisayari]

There are a bug when the latitude is near of 90 or -90 and/or when the longitude is near 180 or -180.

Code to reproduct the bug:

polycircle = polycircles.Polycircle(latitude=0,
                                        longitude=179,
                                        radius=1300*1000,
                                        number_of_vertices=36)
    pol=kml.newpolygon(name="test",outerboundaryis=polycircle.to_kml())
    pol.style.polystyle.color= '40008000' 

EDIT : This bug seems also to be specific to a big radius.

[adamatan]

Thanks. Let me take a look.

[anisayari]

@adamatan I took a look and I think the issue come from geographiclib , what do you think about that ?

[adamatan]

Hi, Apologies for the delayed answer, I haven't had the time to take a deep look yet. I think your assumption is correct. Would you like to write a small code sample that manifests the error and open a bug for geographiclib? Best, and thanks again, Adam

[anisayari]

I run some test directly with geographiclib with this code below and it's work. geograhiclib provide us the right coordinates. So this bug is not from geographiclib

from geographiclib import geodesic

latitude = 30
longitude = 50
radius = 1500 * 1000
number_of_vertices = 30

for i in range(number_of_vertices):
    degree = 360.0 / number_of_vertices * i
    vertex = geodesic.Geodesic.WGS84.Direct(latitude, longitude,
                                                    degree, radius)
    print str(vertex['lon2'])+','+str(vertex['lat2'])

When I use polycircles , the same coordinates are in the kml file with this code :

from polycircles import polycircles
import simplekml

latitude = 30
longitude = 50
radius = 1500 * 1000
number_of_vertices = 30

kml = simplekml.Kml()
polycircle = polycircles.Polycircle(latitude=latitude,longitude=longitude,radius=radius,number_of_vertices=36)
pol=kml.newpolygon(name="test",outerboundaryis=polycircle.to_kml())
pol.style.polystyle.color= '40008000'
kml.save('test.kml')

Output kml :

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
    <Document id="feat_1">
        <Style id="stylesel_0">
            <PolyStyle id="substyle_0">
                <color>40008000</color>
                <colorMode>normal</colorMode>
                <fill>1</fill>
                <outline>1</outline>
            </PolyStyle>
        </Style>
        <Placemark id="feat_2">
            <name>test</name>
            <styleUrl>#stylesel_0</styleUrl>
            <Polygon id="geom_0">
                <outerBoundaryIs>
                    <LinearRing id="geom_1">
                        <coordinates>179.0,43.5166036492,0.0 -177.819217935,43.2745706852,0.0 -174.797978944,42.561338514,0.0 -172.075956287,41.4136364542,0.0 -169.759221823,39.8871192921,0.0 -167.915481841,38.049977521,0.0 -166.577093711,35.9767702286,0.0 -165.748448021,33.7435464003,0.0 -165.414324023,31.4245649224,0.0 -165.547060991,29.0903969861,0.0 -166.111703636,26.8069778508,0.0 -167.069148615,24.6351753897,0.0 -168.377718894,22.6305419376,0.0 -169.993690726,20.8430360534,0.0 -171.871253722,19.3166081319,0.0 -173.962290995,18.0886299048,0.0 -176.21626123,17.1892129739,0.0 -178.580356493,16.6405057005,0.0 179.0,16.4560788387,0.0 176.580356493,16.6405057005,0.0 174.21626123,17.1892129739,0.0 171.962290995,18.0886299048,0.0 169.871253722,19.3166081319,0.0 167.993690726,20.8430360534,0.0 166.377718894,22.6305419376,0.0 165.069148615,24.6351753897,0.0 164.111703636,26.8069778508,0.0 163.547060991,29.0903969861,0.0 163.414324023,31.4245649224,0.0 163.748448021,33.7435464003,0.0 164.577093711,35.9767702286,0.0 165.915481841,38.049977521,0.0 167.759221823,39.8871192921,0.0 170.075956287,41.4136364542,0.0 172.797978944,42.561338514,0.0 175.819217935,43.2745706852,0.0 179.0,43.5166036492,0.0</coordinates>
                    </LinearRing>
                </outerBoundaryIs>
            </Polygon>
        </Placemark>
    </Document>
</kml>

The kml is fine. But when you display it in the kml :

And if you try to change altitude directly from properties , a correct circle is drawing.

The things really weird is that for right coordinates circle, Google earth draw a wrong circle when lat is near of 90/-90 , or and longitude 180/-180. I'm starting to think if it's not an issue directly from Google Earth and kml interpretation, because the coordinates are right.

And if you just try to display this kml without

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2">
    <Document id="feat_1">
        <Style id="stylesel_0">
            <PolyStyle id="substyle_0">
                <color>40008000</color>
                <colorMode>normal</colorMode>
                <fill>1</fill>
                <outline>1</outline>
            </PolyStyle>
        </Style>
        <Placemark id="feat_2">
                    <LineString id="geom_1">
                        <coordinates>179.0,43.5166036492,0.0 -177.819217935,43.2745706852,0.0 -174.797978944,42.561338514,0.0 -172.075956287,41.4136364542,0.0 -169.759221823,39.8871192921,0.0 -167.915481841,38.049977521,0.0 -166.577093711,35.9767702286,0.0 -165.748448021,33.7435464003,0.0 -165.414324023,31.4245649224,0.0 -165.547060991,29.0903969861,0.0 -166.111703636,26.8069778508,0.0 -167.069148615,24.6351753897,0.0 -168.377718894,22.6305419376,0.0 -169.993690726,20.8430360534,0.0 -171.871253722,19.3166081319,0.0 -173.962290995,18.0886299048,0.0 -176.21626123,17.1892129739,0.0 -178.580356493,16.6405057005,0.0 179.0,16.4560788387,0.0 176.580356493,16.6405057005,0.0 174.21626123,17.1892129739,0.0 171.962290995,18.0886299048,0.0 169.871253722,19.3166081319,0.0 167.993690726,20.8430360534,0.0 166.377718894,22.6305419376,0.0 165.069148615,24.6351753897,0.0 164.111703636,26.8069778508,0.0 163.547060991,29.0903969861,0.0 163.414324023,31.4245649224,0.0 163.748448021,33.7435464003,0.0 164.577093711,35.9767702286,0.0 165.915481841,38.049977521,0.0 167.759221823,39.8871192921,0.0 170.075956287,41.4136364542,0.0 172.797978944,42.561338514,0.0 175.819217935,43.2745706852,0.0 179.0,43.5166036492,0.0</coordinates>
                    </LineString>
        </Placemark>
    </Document>
</kml>

You get a good circle.

[LoneWanderer-GH]

Accoding to this stackoverflow link the behavior may differ depending on Google Earth versions. Extrusion to a rather high altitude seems to do the trick. Does anyone have links to some GoogleEarth issue tracker, that could give some clarity to this ?

[adamatan]

We can mention their twitter account. Would you like me to do it? https://twitter.com/googleearth?ref_src=twsrc%5Egoogle%7Ctwcamp%5Eserp%7Ctwgr%5Eauthor

[LoneWanderer-GH]

According to https://developers.google.com/kml/documentation/kmlreference?hl=en we can use longitudes > 180 or < -180 in case there is a anti meridian overlap (at least for overlays).

ˋeastˋ Specifies the longitude of the east edge of the bounding box, in decimal degrees from 0 to ±180. (For overlays that overlap the meridian of 180° longitude, values can extend beyond that range.) `westˋ Specifies the longitude of the west edge of the bounding box, in decimal degrees from 0 to ±180. (For overlays that overlap the meridian of 180° longitude, values can extend beyond that range.)

If I read correctly, the examples above contain values that are within the [-180 ; 180]. It could be interesting to force geographic lib to unroll longitudes (and maybe latitudes) outside of the usual range.

ˋLONG_UNROLL = 32768ˋ Unroll longitudes, rather than reducing them to the range [-180d,180d].

Source https://geographiclib.sourceforge.io/html/python/code.html

[LoneWanderer-GH]

Did a few tests and got it working when crossing antimerdian line. There are still artefacts near poles, but this is related to the latitude (did not found any "LAT_UNROLL" param yet)

The trick is

vertex = geodesic.Geodesic.WGS84.Direct(
    lat1=latitude, lon1=longitude, azi1=degree, s12=radius,
    outmask=geodesic.GeodesicCapability.STANDARD | geodesic.GeodesicCapability.LONG_UNROLL)

KML / screens

test.zip

Code

import itertools
import simplekml
import tqdm
from geographiclib import geodesic
import numpy as np
from tqdm import tqdm
DEFAULT_NUMBER_OF_VERTICES = 36

class Shape(object):
    """Common operations for shapes"""

    def to_lat_lon(self):
        """Returns a tuple of (lat, lon) tuples of the polygon."""
        return self.vertices

    def to_kml(self):
        """Returns a tuple of (lon, lat) tuples suitable for a KML polygon."""
        return self.to_lon_lat()

    def to_wkt(self):
        """Returns a WKT (Well Known Text) representation of the polygon. Note
           that WKT tuples are (lon, lat), not (lat, lon)."""
        vertices = ", ".join(["%f %f" % v for v in self.to_lon_lat()])
        return "POLYGON ((%s))" % vertices

    def to_lon_lat(self):
        """Returns a tuple of (lon, lat) tuples of the polygon.
           The (lon, lat) notation is used in KMLs and WKTs."""
        return tuple(v[::-1] for v in self.vertices)

    def __iter__(self):
        return iter([{'lat': p[0], 'lon': p[1], 'index': self.vertices.index(p)} for p in self.vertices])

    def __str__(self):
        return self.to_wkt()

class Polycircle(Shape):
    """A polygonial approximation of a circle in WGS84 coordinates.
    >>> import polycircles
    >>> number_of_vertices = 20
    >>> polycircle = polycircles.Polycircle(latitude=31.611878, longitude=34.505351, radius=100, number_of_vertices=number_of_vertices)
    >>> len(polycircle.to_lat_lon()) == number_of_vertices + 1
    True
    """

    def __init__(self, latitude, longitude, radius, number_of_vertices=DEFAULT_NUMBER_OF_VERTICES):
        """
        Arguments:
        latitude -- WGS84 latitude, between -90.0 and 90.0.
        longitude -- WGS84 longitude, between -180.0 and 180.0.
        radius -- Circle radius in meters.
        number_of_vertices -- Number of vertices for the approximation
        polygon.
        """
        # Value assertions
        assert number_of_vertices >= 3, "The minimal number of vertices in a polygon is 3."
        assert radius > 0, "Radius can only have positive values."
        assert -180 <= longitude <= 180, "Longitude must be between -180 and 180 degrees."
        assert -90 <= latitude <= 90, "Latitude must be between -90 and 90 degrees."

        self.latitude = latitude
        self.longitude = longitude
        self.radius = radius
        self.number_of_vertices = number_of_vertices

        vertices = []
        for i in range(number_of_vertices):
            degree = 360.0 / number_of_vertices * i
            vertex = geodesic.Geodesic.WGS84.Direct(lat1=latitude, lon1=longitude, azi1=degree, s12=radius,
                                                    outmask=geodesic.GeodesicCapability.STANDARD | geodesic.GeodesicCapability.LONG_UNROLL)
            lat = vertex['lat2']
            lon = vertex['lon2']
            vertices.append((lat, lon))
        vertices.append(vertices[0])
        self.vertices = tuple(vertices)

if __name__ == '__main__':
    # latitude = 45.0
    # longitude = 179.5
    radius = 1500 * 1000
    number_of_vertices = 30
    lat_boundary = 89.89
    lat_step = 10.
    lats = np.arange(-lat_boundary, lat_boundary, lat_step)
    # lon_boundary = 179.89
    # lon_step = 5.0
    # lons = np.arange(-lon_boundary, lon_boundary, lon_step)
    lons = [-178., -179.5, 179.5, 178.]
    # positions = list(itertools.product(lats, lons))
    positions = [(x,y) for x in lats for y in lons]
    kml = simplekml.Kml()
    nb = len(positions)
    for i, pos in tqdm(enumerate(positions), desc="Generating circles for all combinations", total=nb):
        kml_test = kml.newfolder(name=f"Test_{i+1:04_d}/{nb:04_d}")
        center = kml_test.newpoint(name=f"P{i+1}")
        center.coords = [(pos[1], pos[0])]
        polycircle = Polycircle(latitude=pos[0], longitude=pos[1], radius=radius, number_of_vertices=64)
        pol = kml_test.newpolygon(name="test", outerboundaryis=polycircle.to_kml())
        pol.style.polystyle.color = '40008000'
    kml.save('test.kml')

[adamatan]

Fantastic! Would you like to submit a pull request, or do you want me to do it? Adam

On Fri, Dec 3, 2021 at 1:52 PM LoneWanderer-GH @.***> wrote:

Did a few tests and got it working when crossing antimerdian line. There are still artefacts near poles, but this is related to the latitude (did not found any "LAT_UNROLL" param yet)

The trick is

vertex = geodesic.Geodesic.WGS84.Direct( lat1=latitude, lon1=longitude, azi1=degree, s12=radius, outmask=geodesic.GeodesicCapability.STANDARD | geodesic.GeodesicCapability.LONG_UNROLL)

KML / screens

[image: Example_01] https://user-images.githubusercontent.com/19716859/144597748-c1e3bbdc-4507-462f-adea-206e15c5550a.jpg test.zip https://github.com/adamatan/polycircles/files/7649359/test.zip Code

import itertoolsimport simplekmlimport tqdmfrom geographiclib import geodesicimport numpy as npfrom tqdm import tqdmDEFAULT_NUMBER_OF_VERTICES = 36

class Shape(object): """Common operations for shapes"""

def to_lat_lon(self):
    """Returns a tuple of (lat, lon) tuples of the polygon."""
    return self.vertices

def to_kml(self):
    """Returns a tuple of (lon, lat) tuples suitable for a KML polygon."""
    return self.to_lon_lat()

def to_wkt(self):
    """Returns a WKT (Well Known Text) representation of the polygon. Note           that WKT tuples are (lon, lat), not (lat, lon)."""
    vertices = ", ".join(["%f %f" % v for v in self.to_lon_lat()])
    return "POLYGON ((%s))" % vertices

def to_lon_lat(self):
    """Returns a tuple of (lon, lat) tuples of the polygon.           The (lon, lat) notation is used in KMLs and WKTs."""
    return tuple(v[::-1] for v in self.vertices)

def __iter__(self):
    return iter([{'lat': p[0], 'lon': p[1], 'index': self.vertices.index(p)} for p in self.vertices])

def __str__(self):
    return self.to_wkt()

class Polycircle(Shape): """A polygonial approximation of a circle in WGS84 coordinates. >>> import polycircles >>> number_of_vertices = 20 >>> polycircle = polycircles.Polycircle(latitude=31.611878, longitude=34.505351, radius=100, number_of_vertices=number_of_vertices) >>> len(polycircle.to_lat_lon()) == number_of_vertices + 1 True """

def __init__(self, latitude, longitude, radius, number_of_vertices=DEFAULT_NUMBER_OF_VERTICES):
    """        Arguments:        latitude -- WGS84 latitude, between -90.0 and 90.0.        longitude -- WGS84 longitude, between -180.0 and 180.0.        radius -- Circle radius in meters.        number_of_vertices -- Number of vertices for the approximation        polygon.        """
    # Value assertions
    assert number_of_vertices >= 3, "The minimal number of vertices in a polygon is 3."
    assert radius > 0, "Radius can only have positive values."
    assert -180 <= longitude <= 180, "Longitude must be between -180 and 180 degrees."
    assert -90 <= latitude <= 90, "Latitude must be between -90 and 90 degrees."

    self.latitude = latitude
    self.longitude = longitude
    self.radius = radius
    self.number_of_vertices = number_of_vertices

    vertices = []
    for i in range(number_of_vertices):
        degree = 360.0 / number_of_vertices * i
        vertex = geodesic.Geodesic.WGS84.Direct(lat1=latitude, lon1=longitude, azi1=degree, s12=radius,
                                                outmask=geodesic.GeodesicCapability.STANDARD | geodesic.GeodesicCapability.LONG_UNROLL)
        lat = vertex['lat2']
        lon = vertex['lon2']
        vertices.append((lat, lon))
    vertices.append(vertices[0])
    self.vertices = tuple(vertices)

if name == 'main':

latitude = 45.0

# longitude = 179.5
radius = 1500 * 1000
number_of_vertices = 30
lat_boundary = 89.89
lat_step = 10.
lats = np.arange(-lat_boundary, lat_boundary, lat_step)
# lon_boundary = 179.89
# lon_step = 5.0
# lons = np.arange(-lon_boundary, lon_boundary, lon_step)
lons = [-178., -179.5, 179.5, 178.]
# positions = list(itertools.product(lats, lons))
positions = [(x,y) for x in lats for y in lons]
kml = simplekml.Kml()
nb = len(positions)
for i, pos in tqdm(enumerate(positions), desc="Generating circles for all combinations", total=nb):
    kml_test = kml.newfolder(name=f"Test_{i+1:04_d}/{nb:04_d}")
    center = kml_test.newpoint(name=f"P{i+1}")
    center.coords = [(pos[1], pos[0])]
    polycircle = Polycircle(latitude=pos[0], longitude=pos[1], radius=radius, number_of_vertices=64)
    pol = kml_test.newpolygon(name="test", outerboundaryis=polycircle.to_kml())
    pol.style.polystyle.color = '40008000'
kml.save('test.kml')

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/adamatan/polycircles/issues/2#issuecomment-985457369, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHXTLMUZO4H6ECIVZGMTQ3UPCVOBANCNFSM4DFW5PRQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

[LoneWanderer-GH]

You can do it, its a fairly simple fix. Considering the polar region, the issue also exist with the web version of Google Earth. KML Linestrings are ok, only polygons have odd behaviors. When you use the manual polygon tool, the effects are the same when trying to circle the poles.

[adamatan]

Fantastic. I'll try to fix ASAP - a bit busy with a new job right now. Best, Adam

On Fri, Dec 3, 2021 at 5:54 PM LoneWanderer-GH @.***> wrote:

You can do it, its a fairly simple fix. Considering the polar region, the issue also exist with the web version of Google Earth. KML Linestrings are ok, only polygons have odd behaviors. When you use the manual polygon tool, the effects are the same when trying to circle the poles.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/adamatan/polycircles/issues/2#issuecomment-985631475, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHXTLPRL76XZTV24FVXBADUPDR4VANCNFSM4DFW5PRQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.