How can I fit a circle (ideally an ellipse) to latitude/longitude scatter points in Cartopy?

Question

I simply want a best-fitting ellipse or circle for scatter data I have. I have been able to fit a circle to the data using numerous packages, but then the results are clearly nonsense. Maybe I need to do something weird to get results that work (a) for lat/lon data and (b) with Cartopy projections?

I have the following array of longitude/latitude values:

coords = np.array([-153.1906979 ,   62.01707771],
       [  13.05660412,   63.15537447],
       [-175.82610203,   67.11698477],
       [ -10.31730643,   61.74562855],
       [ 168.02402748,   79.60818152],
       [ -34.46162907,   65.10894426],
       [ -57.20962503,   59.49626998],
       [ 113.70202771,   68.22239091],
       [ -80.43411993,   55.6654176 ],
       [  93.77252509,   76.19392633],
       [-104.10892084,   56.68264351],
       [  66.36158188,   67.59664968],
       [-127.75176924,   57.31577071],
       [-151.83057714,   61.64142205],
       [  17.44848859,   56.02194986],
       [-176.30087703,   66.5955554 ],
       [  -5.48747931,   61.95844561],
       [ 160.22917767,   66.07650153],
       [ -27.93440014,   67.82152994],
       [ 137.09393573,   63.71148003],
       [ -53.3290508 ,   55.79699915],
       [ 109.42329666,   75.43090294],
       [ -76.59105583,   59.18143738],
       [  89.94733587,   63.50658353],
       [-100.54585734,   55.16704225],
       [  66.15810397,   64.64851675],
       [-123.65415058,   60.14507524],
       [  41.00262656,   70.67714209],
       [-145.66917977,   68.55315102],
       [  18.34306395,   67.62222778])

I plot them on a map as following:

fig = plt.figure(figsize=(20,20))
ax = fig.add_subplot(121,projection=ccrs.NearsidePerspective(central_longitude=0, central_latitude=90,
                                                             satellite_height=30785831))

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '50m', facecolor='#daf7f7', alpha=0.7, zorder=0))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', facecolor='#ebc7a4', edgecolor='black', alpha=0.7,zorder=0))

ax.set_global()

grid = ax.gridlines(draw_labels=True)
grid.xlabel_style = {'size': 20, 'color': 'black'}
grid.ylabel_style = {'size': 20, 'color': 'black'}


ax.scatter(coords[0:,0], coords[0:,1], c='red', s=40, zorder=1, transform=PlateCarree())

I get this plot

All I now want to do is fit an ellipse or a circle to this scatter data. Using the solution here: https://stackoverflow.com/a/52877062/17583970, I cannot even try plot anything because the b axis of the ellipse is just a nan. Using skg.nsphere_fit() gave a radius of 433, which is obviously wrong or needs transforming in some way.

Any help would be greatly appreciated.

Answer 1

The coordinates you use really matter in this case. Any fitting will use some sort of distance metric that's being minimized, and the distance in lat/lon coordinates doesn't reflect something in meters or miles at all:
https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid

If you're willing to assume that your map projection is a reasonable distance metric you can simply transform the coordinates, and perform the fit on those. I suspect that in this case the fit will be slightly biassed towards the pole (center of the map), making the ellipse a little smaller than it would be when fitted using the actual distance on the earths surface.

Using Cartopy to define the projections of your data and map:

data_proj = ccrs.PlateCarree()
map_proj = ccrs.NearsidePerspective(central_longitude=0, central_latitude=90, satellite_height=30785831)

Those can then be used to convert the coordinates to the map projection:

coords_map = map_proj.transform_points(data_proj, coords[:,0], coords[:,1])[:, :-1]

Fitting (and predicting) an ellipse using Scikit-Image:

from skimage.measure import EllipseModel

model = EllipseModel()
model.estimate(coords_map)
ellipse_coords = model.predict_xy(np.linspace(0, 2*np.pi, 100))

That gives the vertices of the ellipse in the map projection. You could consider converting them back to lat/lon, which would allow you to use ccrs.Geodetic and have Cartopy plot the segments as great circles. But sampling the predicted ellipse in map coordinates might already be fine.

This results in:

fig, ax = plt.subplots(figsize=(6,6), subplot_kw=dict(projection=map_proj), dpi=86, facecolor="w")

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'ocean', '110m', facecolor='#C6EEFF', alpha=1, zorder=0))
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '110m', facecolor='#A2BAA4', edgecolor='none', alpha=1,zorder=0))

ax.set_global()
grid = ax.gridlines(draw_labels=True, lw=0.5, color="k", alpha=0.2)

ax.plot(coords_map[0:,0], coords_map[0:,1], "o", mfc="none", mec='#A70000', mew=1, ms=5, zorder=1, transform=map_proj)
ax.plot(ellipse_coords[:, 0], ellipse_coords[:,1], "-", color="#A70000", transform=map_proj)