Add LatLong "projection"

[dopplershift]

Go ahead and wrap Proj.4's latlong pseudo-projection. This should provide an alternative to PlateCarree to fix #709.

[dopplershift]

Looks non-trivial to make latlong accept something other than radians. Possible with pipelines, but not happening for 0.17.

[deeplycloudy]

Unless there is some other configuration of pyproj being used, it has always been my understanding that latlong takes degrees.

import pyproj as proj4
lon = np.array([-101., -101.])
lat = np.array([33., 34.]) # Should be about 100 km distance
alt= np.array([0.,0.])
ellipse=datum='WGS84'
ERSlla = proj4.Proj(proj='latlong', ellps=ellipse, datum=datum)

# ECEF cartesian, meters
ERSxyz = proj4.Proj(proj='geocent', ellps=ellipse, datum=datum)
x,y,z = proj4.transform(ERSlla, ERSxyz, lon, lat, alt)
dx = np.diff(x)
dy = np.diff(y)
dz = np.diff(z)
print(np.sqrt(dx*dx + dy*dy + dz*dz)/1000.0, ' km')

Which results in 111 km.

[QuLogic]

@dopplershift Do you mean https://proj4.org/operations/conversions/latlon.html? Is that not Geodetic?

>>> ccrs.Geodetic().proj4_init
'+datum=WGS84 +ellps=WGS84 +proj=lonlat +no_defs'

Or would you like a Projection instead of/in addition to a CRS?

[dopplershift]

I want a Projection so that we can hand it to matplotlib. Geodetic does not work for this.

[dopplershift]

So I've made this class just to test:

class LatLong2(_CylindricalProjection):
    def __init__(self, central_longitude=0.0, globe=None):
        proj4_params = [('proj', 'lonlat'), ('datum', 'WGS84')]
        x_max = 180
        y_max = 90
        self._threshold = 0.001

        super(_CylindricalProjection, self).__init__(proj4_params, x_max,
                                                     y_max, globe=globe)

    @property
    def threshold(self):
        return self._threshold

And yet with this code below I need to convert degrees to radians to make it work:

fig = plt.figure()
proj = ccrs.LatLong2()
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines()
ax.plot(-105 * np.pi/180., 40 * np.pi/180., 'ro', transform=proj)
ax.set_extent([-130, -50, 25, 55], crs=ccrs.PlateCarree())

So now that I'm posting this publicly, I await the point at which someone points out how I'm being an idiot. 😀

[deeplycloudy]

Well, I’m the one that spoke without realizing that cartopy has its own proj.4 wrapper, and I was using the pyproj wrapper in my example above, instead of the cartopy _crs.pyx and its class CRS.

So I may be responsible for some noise. But here is some analysis of what pyproj does. It explains why I’ve never worried about radians when using it.

cartopy’s Projection inherits (via metaclass magic) CRS, which does indeed treat the latlong (geodetic) projection with some special deg2rad logic in CRS.transform_point[s]. As the proj.4 API notes point out, pj_transform does want radians.

In pyproj, Proj._transform makes the call to pj_transform, and the arguments include a flag radians used in concert with a check for latlong. The default behavior is radians=False, and this is exposed to the user in pyproj.Proj.transform.

In cartopy, there does not seem to be the option to choose degrees or radians - latlong is always in degrees.

So that is a difference in the proj wrappers. Given my track record in this issue, my comments may be another red herring!

I don’t understand why the example above works the way it does - if matplotlib eventually results in a call to .transform_point[s] and the src_crs is latlong, then it should undergo a deg2rad conversion in CRS.transform_point[s].