Картографическое построение точек неправильно с ортогональной проекцией

Question

Я пытаюсь построить точки на карте, используя Cartopy с Anaconda Python, но получаю некоторые странные ошибки при преобразовании. В моем простом примере я пытаюсь построить 3 очка, но они удваиваются.

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()

Протестировано с:

cartopy==0.16.0 а также Cartopy-0.16.1.dev179-

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

Мой единственный намек был, что с Cartopy-0.16.1.dev179- а также proj4==5.0.1Я получил это 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.

Я открыл проблему на https://github.com/SciTools/cartopy/issues/1172 но проблема была закрыта. Кто-нибудь знает, как заставить картопи работать правильно с ортогональными проекциями?

Answer 1

Насколько я знаю, есть несколько подходов, которые вы могли бы использовать, чтобы получить ожидаемый результат.

Во-первых, явным образом преобразуйте точки, чтобы быть в родной проекции...

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()

Это сюжеты, как и ожидалось...

Ожидаемый орфографический проекционный участок

Оригинальная проблема, которую вы видите, в том, что cartopy пытается интерпретировать последовательность точек как ограниченную геометрию (или путь), но становится немного запутанным. Явное преобразование точек широты / долготы в родные точки орфографии уклоняется от этой пули.

Зная этот кусок информации, мы могли бы альтернативно вызвать соответствующий метод, который учитывает список точек как отдельные точки (и избегать cartopy делать предположения, которые не соответствуют нашим ожиданиям) с помощью scatter вместо 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()

Это также работает для меня.

НТН