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

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

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

НТН