How do you display the scale in meters, the north arrow and the axes in latitude and longitude on a map with Geopandas?

With reference to this issue, is it possible to have the scale bar (projected in meters, so 3857 for example) with the x,y axes in latitude, longitude projection (4326) and the north arrow?

I don't see a turnkey solution to do this with geopandas. While this seems to be basic settings for map display with GIS. Is there a technical reason for this?

import geopandas as gpd
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt

df = gpd.read_file(gpd.datasets.get_path('nybb'))
ax = df.to_crs(4326).plot()
ax.add_artist(ScaleBar(1)) #how add ScaleBar for df in 3857? 
plt.show()

From this, it looks like you have to compute the great circle distance between two locations A and B with coordinates A=[longitudeA,latitudeA] and B=[longitudeA+1,latitudeA], at the latitude you are interested in (in your case ~40.7°). To compute the great circle distance you can use the 'haversine_distances' from sklearn (here) and multiply it by the radius of the earth 6371000 to get the distance in meters. Once you get this distance dx, you can just pass it to your scalebar with ScaleBar(dx=dx,units="m").

So overall, the code looks like that:

import numpy as np
import geopandas as gpd
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import haversine_distances

df = gpd.read_file(gpd.datasets.get_path('nybb'))
ax = df.to_crs(4326).plot()
A=[-74.5*np.pi/180.,40.7*np.pi/180.] #Latitude of interest here 40.7 deg, longitude -74.5
B=[-73.5*np.pi/180.,40.7*np.pi/180.] ##Latitude of interest here 40.7 deg, longitude -74.5+1
dx=(6371000)*haversine_distances([A,B])[0,1]
ax.add_artist(ScaleBar(dx=dx,units="m")) 
plt.show()

And the output gives: