Cartopy plotting points incorrectly with Orthographic projection

user_1945087 0
Viewed 2602

I am trying to plot map points using Cartopy with Anaconda Python but am getting some strange failures with the transform. In my simple example, I am trying to plot 3 points but they are getting doubled.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

lons = [214.5, 2.7, 197.5]                                                                                               
lats = [35, 36, 37.]

ax = plt.axes(projection=ccrs.Orthographic(
    central_longitude=0,
central_latitude=90))
# plot lat/lon points                                                                                                         
ax.plot(lons, lats, 'ro',
        transform=ccrs.Geodetic())
# plot north pole for reference                                                                                               
ax.plot([0], [90], 'b^',
    transform=ccrs.Geodetic())
# add coastlines for reference                                                                                                
ax.coastlines(resolution='50m')
ax.set_global()

plt.show()

Plot output

Tested with:

cartopy==0.16.0 and Cartopy-0.16.1.dev179-

proj4==4.9.3, proj4==5.0.1, proj4==5.0.2

My only hint was that with Cartopy-0.16.1.dev179- and proj4==5.0.1, I got this UserWarning:

/Users/***/anaconda3/lib/python3.6/site-packages/cartopy/crs.py:1476: UserWarning: The Orthographic projection in Proj between 5.0.0 and 5.1.0 incorrectly transforms points. Use this projection with caution.

I opened an issue on https://github.com/SciTools/cartopy/issues/1172 but the issue was closed. Anyone know how to get cartopy working correctly with Orthographic projections?

cartopy python
1 Answers

As far as I know, there are a couple of approaches that you could use to get the result that you expect.

Firstly, explicitly transform the points to be in the native projection...

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# create the lat/lon points
lons = np.array([214.5, 2.7, 197.5])
lats = np.array([35, 36, 37.])

# create the projections
ortho = ccrs.Orthographic(central_longitude=0, central_latitude=90)
geo = ccrs.Geodetic()

# create the geoaxes for an orthographic projection
ax = plt.axes(projection=ortho)

# transform lat/lons points to othographic points
points = ortho.transform_points(geo, lons, lats)

# plot native orthographic points                                                                                
ax.plot(points[:, 0], points[:, 1], 'ro')

# plot north pole for reference (with a projection transform)                                                                                           
ax.plot([0], [90], 'b^', transform=geo)

# add coastlines for reference                                                                                                
ax.coastlines(resolution='50m')
ax.set_global()

This plots as expected...

Expected Orthographic projection plot

The original issue that you're seeing is that cartopy is attempting interpret the sequence of points as a bounded geometry (or path), but is getting a little confused. Explicitly converting the lat/lon points to be native orthographic points dodges this bullet.

Knowing this nugget of information, we could alternatively call the appropriate method that respects the list of points as individual points (and avoid cartopy making assumptions that don't meet our expectations) by using scatter instead of plot...

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# create the lat/lon points
lons = np.array([214.5, 2.7, 197.5])
lats = np.array([35, 36, 37.])

# create the projections
ortho = ccrs.Orthographic(central_longitude=0, central_latitude=90)
geo = ccrs.Geodetic()

# create the geoaxes for an orthographic projection
ax = plt.axes(projection=ortho)

# plot native orthographic scatter points                                                                                
ax.scatter(lons, lats, marker='o', c='r', transform=geo)

# plot north pole for reference                                                                                               
ax.plot([0], [90], 'b^', transform=geo)

# add coastlines for reference                                                                                                
ax.coastlines(resolution='50m')
ax.set_global()

This also works for me.

HTH

user_17577870
Related