Allow user to control transform resolution

[pelson]

Currently, the great circle transform resolution has been set to an arbitrary value. Tweak the value to some heuristic, and make it possible for users to define the value themselves.

[pelson]

@rhattersley : the threshold property on a CRS is presumably designed to do this? What are your intentions for that?

A good example of the problems that are currently being seen:

import matplotlib.pyplot as plt
import numpy

import cartopy.crs as ccrs


prj = ccrs.RotatedPole(180, 37)
ax = plt.axes(projection=prj)
ax.coastlines()
ax.gridlines()

v = 10
ax.set_extent([-v-5, v+5, 45, 60])
plt.plot([-v, v], [50, 50], transform=ccrs.PlateCarree(), color='blue', lw=2, zorder=10)
plt.plot([-v, v], [50, 50], transform=ccrs.Geodetic(), color='red', lw=5)

plt.show()

Tweaking the value of v shows that the threshold, as it stands, does not interpolate suitably for limited area maps. (the geodetic line is supposed to be fairly straight due to the pole rotation).

[rhattersley]

The original tolerance values were set back when the interpolation was a lot slower, so were heavily weighted to bias performance over accuracy. So they could probably be set tighter now. (For the display of great circles that might well be enough.)

But that said, I'd always hoped we could make the tolerance values "context aware", so they would adjust to suit the current plot. Simplifying somewhat, we want to put in enough points on the interpolation to ensure the interpolated curve is accurate to the order of a pixel - no more, no less.

[pelson]

FWIW an ugly workaround is to threshold hack the map projection:

class RP(ccrs.RotatedPole):
    @property
    def threshold(self):
        return 0.01


prj = RP(180, 37)

In general, we have a few read-only thresholds that we should really make easy to override until we come up with a better architecture to adaptively set the threshold.

[Guymer]

I would love for the transform resolution to become a user-accessible value too. I completely understand, and support, @rhattersley 's comments about the reasons for it being hard-coded ages ago. However, I too am finding that large great circles (when viewed on a limited area map) are very jagged.

[aurelgriesser]

This was opened in 2012, is anyone still interested in this?
I just used @pelson's suggested (5 Sep 2017) custom class to make an AlbersEqualArea map correctly draw a great circle line. Seems to work but would be nice to have this resolved properly.

[dopplershift]

I'd wager people are still interested, but no one so far has been interested enough to contribute a PR implementing a writable threshold value and making sure that resets everything necessary.

[jonnyhtw]

defo still interested :)

[blazing216]

Still interested! I traced back to this issue when trying to plot a smooth great circle path using cartopy. I took some time to explore the idea of 'context aware threshold' of @rhattersley, seems like it may be easy to implement now. Here is the result of my experiment (using @pelson's example):

import matplotlib.pyplot as plt
import numpy as np

import cartopy.crs as ccrs

dpi = 100


v = 80
def plot_example(threshold):
    ax.coastlines()
    ax.gridlines()
    ###########################################
    # the threshold is changed before plotting the two lines
    ###########################################
    prj.threshold = threshold
    ax.set_extent([-v-5, v+5, 5, 60])
    plt.plot([-v, v], [50, 50], transform=ccrs.PlateCarree(), color='blue', lw=1)
    plt.plot([-v, v], [50, 50], transform=ccrs.Geodetic(), color='red', 
        lw=1)

# I added @property.setter to ccrs.RotatedPole to set threshold
prj = ccrs.RotatedPole(180, 37)


plt.figure(dpi=dpi)

# threshold = 0.5
ax = plt.subplot(211, projection=prj)
plot_example(threshold=0.5)
plt.title('Threshold = 0.5')

# adaptive threshold
ax = plt.subplot(212, projection=prj)

# compute the length of one pixel in the (projected) data coordinates
pixel2data = ax.transData.inverted()
xy0 = pixel2data.transform((0,0))
xy1 = pixel2data.transform((1,1))
#dxy = np.sqrt((xy0[0] - xy1[0])**2 + (xy0[1] - xy1[1])**2)
pixelsize_in_data_coordinates = np.abs(xy0[0] - xy1[0])

plot_example(threshold=pixelsize_in_data_coordinates)
plt.title(f'Adaptive Threshold ({pixelsize_in_data_coordinates:.3g})')

plt.suptitle(f'{dpi} dpi')

plt.savefig(f'test_gcp_clean_{dpi}_dpi.png', dpi=dpi)
#plt.show()