Validity of LinearRing in shapely

Question

Consider the following code:

import numpy as np
from shapely.geometry.polygon import LinearRing
llpts = np.asarray([[18.5165    , 73.82883333],
       [18.51678333, 73.8279    ],
       [18.5156    , 73.82848333],
       [18.51586667, 73.82755   ],
       [18.5165    , 73.82883333]])
linring_ll = LinearRing(llpts)
print('linring_ll.is_valid =',linring_ll.is_valid)

Using Shapely 1.7.1 and Python 3.9.2 the above code prints 'True' on my computer even though it seems that the LinearRing crosses itself (I plotted it in Spyder). Is this a rounding / numerics issue? This goes away when I convert the points from lat/lon space to xy space using equirectangular projection.

The complete code follows:

import numpy as np
from shapely.geometry.polygon import LinearRing, LineString, Polygon
def convert_llpts_xypts(latlon,rad_earth,center_lat,projtype='equirectangular'):
    # convert from latlon space to xy space 
    # latlon must be in degrees
    # lat lon looks like [(lat1,lon1),(lat2,lon2)....]
    # lat lon needs to be a np array
if (projtype == 'equirectangular'):
    # https://stackoverflow.com/questions/16266809/convert-from-latitude-longitude-to-x-y
    _radconv          = np.pi/180.0  
    _rade_radconv     = rad_earth*_radconv
    _rade_radconv_cos = _rade_radconv*np.cos(center_lat*_radconv)

    _xx = _rade_radconv_cos*latlon[:,1]
    _yy = _rade_radconv*latlon[:,0]
    _xx = _xx[:,None]
    _yy = _yy[:,None]

    return np.hstack((_xx,_yy))

else:
    raise ValueError(f'Invalid projtype={projtype} in convert_llpts_xypts')


points in lat/lon space
llpts = np.asarray([[18.5165    , 73.82883333],
       [18.51678333, 73.8279    ],
       [18.5156    , 73.82848333],
       [18.51586667, 73.82755   ],
       [18.5165    , 73.82883333]])
compute a latitude close to the center of the polygon
center_lat = np.sum(llpts[0:4,0])/4.0
convert points to xy-space
xypts = convert_llpts_xypts(llpts,6378137.0,center_lat)
make LinearRings
linring_ll = LinearRing(llpts)
linring_xy = LinearRing(xypts)
print(f'validity of ring in lat/lon space = ', linring_ll.is_valid)
print(f'validity of ring in xy space = ', linring_xy.is_valid)

Answer 1

This is a bug in the underlying GEOS library, and its upstream parent JTS. Reported as JTS-736