Inconsistent behavior when supplying globes to PlateCarree

[samwisehawkins]

Supplying a spherical globe, or even a default ccrs.Globe() to PlateCarree changes the behavior compared to not supplying any globe. I've been transforming latitude and longitude into an LCC projection from a PlateCaree with (a) the WGS84 ellipsoid and (b) a spherical Earth. The results differ by several orders of magnitude. Perhaps this is related to #653?

import cartopy
import cartopy.crs as ccrs

lon,lat = np.meshgrid(np.arange(-5,5,1),np.arange(50,60,1))
sphere = ccrs.Globe(semimajor_axis=6370000, semiminor_axis=6370000, ellipse=None) # have     also tried 'sphere'
ellipsoid = ccrs.Globe()
lcc = ccrs.LambertConformal(central_longitude=0, 
                        central_latitude=55.0, 
                        standard_parallels=(53,57)) 

pc_wgs1= ccrs.PlateCarree()
pc_wgs2= ccrs.PlateCarree(globe=ellipsoid)
pc_sphere = ccrs.PlateCarree(globe=sphere)

r1 = lcc.transform_points(pc_wgs1, lon, lat)
r2 = lcc.transform_points(pc_wgs2, lon, lat)
r3 = lcc.transform_points(pc_sphere, lon, lat)

print(r1[0,0,0])
print(r2[0,0,0])
print(r3[0,0,0])

Gives:

-359276.021664
-7.36249036419
-7.37189510538

[pelson]

Hi @samwisehawkins, thanks for taking the time to give a good example of the problem.

First things first - the inputs for transform_points with a PlateCarree projection are emphatically not latitudes and longitudes - they are the Cartesian coordinates which represent the PlateCarree projection. If we choose a suitable ellipsoid for PlateCarree, then we can make it look very much like the inputs are latitudes and longitudes (which we have decided to do in cartopy). But we should keep in mind that they really are not latitudes and longitudes (this is the reason we get the dotted straight line in http://scitools.org.uk/cartopy/docs/latest/matplotlib/intro.html#adding-data-to-the-map).

If we want a coordinate system which does speak latitudes and longitudes (always) we can look at the Geodetic CRS:

>>> g0 = ccrs.Geodetic()
>>> g1 = ccrs.Geodetic(globe=ellipsoid)

>>> print(lcc.transform_points(g0, lon, lat)[0, 0, 0])
-359276.021664
>>> print(lcc.transform_points(g1, lon, lat)[0, 0, 0])
-359276.021664

There is definitely room for improvement in the cartopy docs on this subject. Any ideas where you would put them @samwisehawkins?

[samwisehawkins]

Hi @pelson, thanks for the info.

I agree I was using PlateCarree rather sloppily there - a widespread problem where it is often equated as being "unprojected". At the risk of opening up a Pandora's box, I think it would be nice to have some general documentation about coordinate systems, datums and projections, probably before "Coordinate reference systems in Cartopy". Perhaps Cartopy is not the place for this kind of background information, but in my experience much documentation found on the internet is imprecise and often misleading.

I'm happy to contribute to those docs.

[dopplershift]

@pelson Having just lost a couple hours to this--because I forgot about passing globe being an issue--I think this is a design mistake, having such widely different behavior when passing in the Globe. No other projection in the library changes so much as to break everything when you change from the default globe, at least in terms of making previously fine coordinate values completely non-sensical. I think two separate concepts are being conflated in the current implementation, and they should be separated:

1. Deprecate/remove the ability to pass a custom globe to PlateCarree
2. Add a EquidistantCylindrical (or Equirectanguler) class that defaults to WGS84.

At that point PlateCarree might even be a simple subclass of EquidistantCylindrical. Thoughts?

edit: I see now this was kind of proposed in #653, so I think that lends support to the idea. To be clear, I am volunteering to make the PR if it's agreed that this is a good idea.