Масштаб проекции картопии не соответствует

Я использую cartopy для отображения KDE, наложенного на карту мира. Первоначально я использовал проекцию ccrs.PlateCarree без проблем, но в тот момент, когда я попытался использовать другую проекцию, мне показалось, что масштаб проекции резко увеличился. Для справки я привел пример, который вы можете протестировать на своем компьютере ниже (просто закомментируйте две строки projec, чтобы переключаться между проекциями)

from scipy.stats import gaussian_kde
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

projec = ccrs.PlateCarree()
#projec = ccrs.InterruptedGoodeHomolosine()


fig = plt.figure(figsize=(12, 12))



ax = fig.add_subplot(projection=projec)

np.random.seed(1)

discrete_points = np.random.randint(0,10,size=(2,400))

kde = gaussian_kde(discrete_points)
x, y = discrete_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
resolution = 1
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]
earth = plt.cm.gist_earth_r

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', 
                                                edgecolor='black', facecolor="none"))

ax.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,
    alpha=0.8,
    cmap=earth, transform=projec)

ax.axis('on')
ax.get_xaxis().set_visible(True)
ax.get_yaxis().set_visible(True)

ax.set_xlim(-30, 90)
ax.set_ylim(-60, 60)

plt.show()

Вы заметите, что при использовании проекции ccrs.PlateCarree() KDE удобно размещается над Африкой, однако при использовании проекции ccrs.InterruptedGoodeHomolosine() вы вообще не видите карту мира. Это потому, что карта мира имеет огромный масштаб. Ниже приведены изображения обоих примеров:

Проекция пластины Карри:  Тарелка Карри проекция

Прерванная проекция Гуда Гомолосин (стандартное увеличение):  без увеличения

Прерванная проекция Гуда Гомолосина (уменьшено):

zoom goode

Если бы кто-нибудь мог объяснить, почему это происходит, и как это исправить, чтобы я мог нанести одни и те же данные на разные проекции, это был бы очень признателен.

РЕДАКТИРОВАТЬ:

Я также хотел бы указать, что я пытался добавить transform=projec в строку 37 в примере, который я включил, а именно:

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', 
                                                edgecolor='black', facecolor="none", transform=projec))

Однако это не помогло. Фактически, казалось, что после добавления этого мира карта мира больше не появлялась.

РЕДАКТИРОВАТЬ:

В ответ на ответ JohanC я получаю следующий сюжет при использовании этого кода:

Сюжет JohanC увеличен

И уменьшено:

Сюжет JohanC уменьшен


person Recessive    schedule 13.11.2020    source источник


Ответы (1)


Комментарии к вашим участкам:

Участок1: (справочная карта)

  • проекция: PlateCarree projection
  • Размер изображения (Zgrid) покрывает (приблизительную) квадратную область, примерно 40 градусов с каждой стороны
  • нижний левый угол изображения находится в широте / долготе: (0,0)

Сюжет 2

В: Почему топографические объекты не отображаются на карте?

О: Сюжет занимает очень маленькую территорию, на которой нет ни одной из них.

  • проекция: Прерванный
  • данные изображения, Zgrid объявляется соответствующим координатам сетки (проекции) (единица измерения: метры)
  • карта строится в пределах нескольких метров как по x, так и по y, и соотношение сторон не равно.

Сюжет 3

В: Почему изображение Zgrid не отображается на карте?

О: График покрывает очень большую область, поэтому изображение становится слишком маленьким для печати.

  • проекция: Прерванная проекция ГудГомолозина
  • экстент изображения (Zgrid) очень маленький, не виден в этом масштабе
  • карта строится с большим размером, и соотношение сторон не одинаковое.

Средства правовой защиты (для Участков 2 и 3)

  • Zgrid требует правильного преобразования из широты / долготы в координаты проекции осей
  • экстенты карты также необходимо преобразовать и установить соответствующим образом
  • соотношение сторон должно быть установлено равным, чтобы предотвратить неравномерное растяжение по осям x и y.

О графиках с линиями сетки

  • полезно для справки о местоположении
  • широта / параллели: ОК с InterruptedGoodeHomolosine в этом случае
  • долгота / меридианы: проблематично (не знаю, как исправить !!)

Вот модифицированный код, который запускается и создает требуемую карту.

# proposed code

from scipy.stats import gaussian_kde
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fig = plt.figure(figsize=(7, 12))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())

np.random.seed(1)
discrete_points = np.random.randint(0,10,size=(2,400))

kde = gaussian_kde(discrete_points)
x, y = discrete_points
# https://www.oreilly.com/library/view/python-data-science/9781491912126/ch04.html
resolution = 1
x_step = int((max(x)-min(x))/resolution)
y_step = int((max(y)-min(y))/resolution)
xgrid = np.linspace(min(x), max(x), x_step+1)
ygrid = np.linspace(min(y), max(y), y_step+1)
Xgrid, Ygrid = np.meshgrid(xgrid, ygrid)
Z = kde.evaluate(np.vstack([Xgrid.ravel(), Ygrid.ravel()]))
Zgrid = Z.reshape(Xgrid.shape)

ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]
earth = plt.cm.gist_earth_r


ocean110 = cfeature.NaturalEarthFeature('physical', 'ocean', \
        scale='110m', edgecolor='none', facecolor=cfeature.COLORS['water'])
ax.add_feature(ocean110, zorder=-5)

land110 = cfeature.NaturalEarthFeature('physical', 'land', '110m', \
        edgecolor='black', facecolor="silver")
ax.add_feature(land110, zorder=5)

# extents used by both Zgrid and axes
ext = [min(x)*5, max(x)*5, min(y)*5, max(y)*5]

# plot the image's data array
# note the options: `extent` and `transform`
ax.imshow(Zgrid,
    origin='lower', aspect='auto',
    extent=ext,  #set image's extent
    alpha=0.75,
    cmap=earth, transform=ccrs.PlateCarree(),
    zorder=10)

# set the plot's extent with proper coord transformation
ax.set_extent(ext, ccrs.PlateCarree())

ax.coastlines()
#ax.add_feature(cfeature.BORDERS) #uncomment if you need

ax.gridlines(linestyle=':', linewidth=1, draw_labels=True, dms=True, zorder=30, color='k')
ax.set_aspect('equal')  #make sure the aspect ratio is 1

plt.show()

Выходная карта:

interup_proj

person swatchai    schedule 14.11.2020
comment
Спасибо за подробный ответ! По большей части это работает, единственная проблема, с которой я столкнулся, - это то, что KDE примерно в 10 раз больше, чем должен быть. В свою правку я включил пример. Я напрямую копирую и вставляю ваш код в сторону от строки ax.gridlines, потому что с моей версией картографии я не могу рисовать их ни на одной другой проекции, кроме PlateCarree или Mercator. - person Recessive; 16.11.2020