Как выглядит рельеф Марса? Выясняем с помощью Python

К старту флагманского курса по Data Science делимся визуализациями марсианского ландшафта на основе изображений, полученных благодаря беспилотным полётам над поверхностью планеты. За подробностями приглашаем под кат.

Лучший способ узнавать что-то новое — решать интересные или необычные задачи. Одна из них — создание топографических карт Марса с помощью сложных цветовых карт.

Благодаря шумихе, поднятой одним американским миллиардером, Марс всё чаще появляется в новостях, и наверняка уже скоро люди будут ходить по его поверхности. Многие представляют Марс большим красным шаром в небе, но мало кто знает, как он выглядит на самом деле. А ведь у Марса богатый ландшафт — от каньонов и кратеров до вулканов.

В ходе многочисленных беспилотных полётов карту детализировали. И теперь легко получить реальное представление о том, как выглядит Марс. Создадим карту планеты и с помощью магии цветовых карт представим, как выглядел бы терраформированный Марс, его океаны и реки. А главное — поймём, как правильно использовать цветовые карты для решения конкретных задач.

Изучение данных

Сначала загружаем данные о ландшафте Марса. Данные правительства США доступны здесь и не имеют никаких ограничений для пользователей. Имейте в виду: объём данных огромен — это 2,3 Гб.

Формат данных — tif, загружаются они через rasterio в np.array. Размеры массива: 23040 x 46090. Диапазон значений — это высоты Красной планеты. На Марсе нет океанов, поэтому нет точки отсчёта высоты — уровня моря. Если коротко, «уровень моря» на Марсе определяется как средняя высота на планете. Отрицательные значения рельефа — это области ниже средней высоты, а положительные — выше.

Гора Олимп — самая высокая на планете (21 241 м относительно средней высоты Марса), а равнина Эллада — самая низкая точка — это –8201 м относительно средней высоты Марса. Разница между самой высокой и самой низкой точками на Марсе (~30 км) больше, чем между Эверестом и Марианской впадиной на Земле (~20 км):

      import rasterio
      import numpy as np

      mars = rasterio.open('../../Planets/mars/data/Mars_MGS_MOLA_DEM_mosaic_global_463m.tif')
      mars = mars.read()
      
      print(mars.shape)
      print(np.amin(mars[0]))
      print(np.amax(mars[0]))
      print(np.amax(mars[0]) + abs(np.amin(mars[0])))

Вывод кода:

      (1, 23040, 46080)
      -8201
      21241
      29442

Нанесённые на карту данные дают представление о том, с чем мы имеем дело. Предупреждение: на стандартном ноутбуке выполнение займёт немало времени. Можно успеть выпить чаю.

      import matplotlib.pyplot as plt

      fig, ax = plt.subplots()
      fig.set_size_inches(14, 7)

      ax.imshow(mars[0], cmap="viridis")
      ax.axis('off')

      plt.show()

При использовании стандартной цветовой карты возникает проблема: отображая на ней значения рельефа линейно, получаются размытые объекты и большинство из них теряется. На Марсе объекты в основном расположены между –8000 и 12 000 м, а средняя высота намного ближе к равнине Элладе, чем к горе Олимп. 

Четыре горы ближе к жёлтому и занимают крошечную по площади часть планеты, но огромны по высоте, и у них четверть всех цветов на цветовой карте. Остальные объекты отображаются в небольшом диапазоне цветов, чёткость плохая:

Определение цветовых карт

Решим проблему, создав цветовую карту и нарушив правила визуализации данных. Нормализуем линейную цветовую карту между наименьшими и наибольшими значениями данных. Но в этом случае всего четыре объекта выше 14 000 м и куча объектов ниже. Поэтому создадим цветовую карту для отображения цветов между –8201 (самая низкая точка в данных) и 14 000 м, а всё, что выше, будет в последнем цвете цветовой карты. matplotlib не особо хороша в отображении рельефа на цветовых картах. Вот что для этого используют в NASA, и мы сделаем то же самое.

Ниже показан раздел кода, в котором создаётся подобная цветовая карта, нормализованная между –8210 и 14 000 м в 10-метровых полосках. Метод LinearSegmentedColormap принимает список цветов, создаётся цветовая карта с N цветами, линейно расположенными между цветами из списка. 

Класс BoundaryNorm позволяет сопоставлять на ней цвета с набором заданных значений. В нашем случае создаётся объект BoundaryNorm, отображающий цвета в 10-метровых полосках между –8210 и 14 000 м. Например, значения от 0 до 10 м сопоставляются с одним цветом, значения от 10 до 20 м — с другим и т. д. Всё, что выше 14 000 м, отобразится в последнем цвете карты:

from matplotlib.colors import BoundaryNorm, LinearSegmentedColormap
import numpy as np

custom_cmap = LinearSegmentedColormap.from_list('mars', ['#162252',
                                                         '#104E8B',
                                                         '#00B2EE',
                                                         '#00FF00',
                                                         '#FFFF00',
                                                         '#FFA500',
                                                         '#FF0000',
                                                         '#8b0000',
                                                         '#964B00',
                                                         '#808080',
                                                         '#FFFFFF'], N=2221)

bounds = np.arange(-8210, 14000, 10)
norm = BoundaryNorm(bounds, custom_cmap.N)

Получилась цветовая карта с примерными высотами марсианских объектов. Цветов стало гораздо больше, ведь диапазон значений увеличился:

Составим карту марсианского рельефа с нормализованными значениями и новой цветовой картой. Просто удивительно: это совсем другая карта, а ведь мы изменили только цветовую карту и способ сопоставления цветов со значениями:

fig, ax = plt.subplots()
fig.set_size_inches(14, 7)

ax.imshow(mars[0], cmap=custom_cmap, norm=norm)
ax.axis('off')
logo = plt.imread('../../Branding/globe_black.png')
newax = fig.add_axes([0.82, 0.13, 0.08, 0.08], anchor='NE')
newax.imshow(logo)
newax.axis('off')
txt = ax.text(0.0, 0.02, "Martian Topography \n@PythonMaps",
              size=4,
              color='white',
              transform = ax.transAxes)
plt.show()

Добавим обозначения всех объектов:

def label_features(ax):

     ax.text(5000, 8400, "Olympus\n Mons")
     ax.text(7000, 10000, "Ascreaus\n Mons")
     ax.text(6200, 11500, "Pavonis\n Mons")
     ax.text(5000, 13000, "Arsia\n Mons")
     ax.text(8500, 13000, "Tharsis", rotation=45)
     ax.text(12000, 14000, "Vallies\n Marineris")
     ax.text(17000, 18000, "Argyre")
     ax.text(20000, 2000, "Vastitas Borealis")
     ax.text(17000, 8500, "Chryse\n Planitia")
     ax.text(20000, 5000, "Acidalia Planitia")
     ax.text(30000, 16500, "Hellas", color='white')
     ax.text(37000, 5000, "Utopia\n Planitia")
     ax.text(41000, 8000, "Elysium")
     ax.text(34000, 9800, "Isidis")
     ax.text(7000, 5000, "Alba Patera")
     ax.text(2000, 6500, "Amazonis\n Planitia")

     return ax

fig, ax = plt.subplots()
fig.set_size_inches(14, 7)

ax.imshow(mars[0], cmap=custom_cmap, norm=norm)
ax = label_features(ax)
ax.axis('off')
newax = fig.add_axes([0.82, 0.13, 0.08, 0.08], anchor='NE')
newax.imshow(logo)
newax.axis('off')
txt = ax.text(0.0, 0.02, "Martian Topography \n@PythonMaps",
              size=4,
              color='white',
              transform = ax.transAxes)
plt.show()

Терраформирование Марса

На самом деле терраформировать Марс мы не будем, а создадим цветовые карты, подобные земным. И представим, каким будет полностью терраформированный Марс. Это цветовые карты с тёмно-синими и голубыми цветами для значений ниже уровня моря и зелёно-коричневыми — для значений выше уровня моря.

В Matplotlib для этого есть цветовая карта terrain, но её синяя часть недостаточно тёмная. Есть ещё карта ocean. Объединив их, мы получим цветовую карту, подобную земным картам. Цель объединения — создать объект BoundaryNorm, чтобы сопоставлять значения ниже уровня моря с картой ocean, а значения выше уровня моря — с картой terrain.

При создании цветовой карты цвета индексируются линейно по шкале от 0 до 1. С помощью списка индексов из цветовой карты извлекается список цветов и исключается её часть. Например, ниже в plt.cm.ocean(np.linspace(0.2, 0.8, 821)) из цветовой карты ocean возвращается список с 821 цветом и исключаются цвета из её части ниже 0,2 и выше 0,8. 

То же самое сделано для цветовой карты terrain, но здесь возвращается 1400 цветов. Затем эти два списка объединяются в одну цветовую карту с помощью метода LinearSegmentedColourmap.

Рельеф будет отображаться на цветовой карте 821 цветом ocean и 1400 цветами terrain. Класс BoundaryNorm позволяет сопоставлять цвета карты с набором заданных значений. 

В коде ниже создадим цветовую карту с 2221 цветом (821 ocean + 1400 terrain) и объектом BoundaryNorm, с помощью которого сопоставим значения данных с цветовой картой по предопределённым уровням — 821 уровень ниже 0 и 1400 выше 0. Все значения данных ниже 0 будут в цветах части цветовой карты ocean, а все значения выше нуля — в цветах карты terrain:

colors_undersea = plt.cm.ocean(np.linspace(0.2, 0.8, 821))
undersea_map = LinearSegmentedColormap.from_list('undersea_map', colors_undersea, N=821)

colors_land = plt.cm.terrain(np.linspace(0.25, 1, 1400))
land_map = LinearSegmentedColormap.from_list('land_map', colors_land, N=1400)

colors = np.vstack((colors_undersea, colors_land))
terrain_map = LinearSegmentedColormap.from_list('cut_terrain', colors, N=2221)

bounds = np.arange(-8210, 14000, 10)
norm = BoundaryNorm(bounds, terrain_map.N)

Данные этой новой цветовой карты с границами цветов делают Марс очень похожим на Землю:

def label_features_terrain(ax):
    ax.text(5000, 8400, "Olympus\n Mons")
    ax.text(7000, 10000, "Ascreaus\n Mons")
    ax.text(6200, 11500, "Pavonis\n Mons")
    ax.text(5000, 13000, "Arsia\n Mons")
    ax.text(8500, 13000, "Tharsis", rotation=45)
    ax.text(12000, 14000, "Vallies\n Marineris")
    ax.text(17000, 18000, "Argyre", color='white')
    ax.text(20000, 2000, "Vastitas Borealis", color='white')
    ax.text(17000, 8500, "Chryse\n Planitia", color='white')
    ax.text(20000, 5000, "Acidalia Planitia", color='white')
    ax.text(30000, 16500, "Hellas", color='white')
    ax.text(37000, 5000, "Utopia\n Planitia", color='white')
    ax.text(41000, 8000, "Elysium")
    ax.text(34000, 9800, "Isidis", color='white')
    ax.text(7000, 5000, "Alba Patera")
    ax.text(2000, 6500, "Amazonis\n Planitia", color='white')
    return ax

fig, ax = plt.subplots()
fig.set_size_inches(14, 7)

ax.imshow(mars[0], cmap=terrain_map, norm=norm)
ax = label_features_terrain(ax)
ax.axis('off')
newax = fig.add_axes([0.82, 0.13, 0.08, 0.08], anchor='NE')
newax.imshow(logo)
newax.axis('off')
txt = ax.text(0.0, 0.02, "Martian Topography \n@PythonMaps",
              size=4,
              color='white',
              transform = ax.transAxes)
plt.show()

Изменение уровня моря

В качестве уровня моря мы выбрали произвольное значение — среднюю высоту Марса. Но этот уровень можно менять, применяя описанный выше метод создания цветовых карт. Цветовая карта создаётся из 2221 цвета: 821 из карты ocean и 1400 из карты terrain. Меняя количество цветов в каждой карте, но сохраняя их общее число, будем менять уровень моря. Там, где на карте заканчивается синий и начинается зелёный.

В функции ниже принимается значение уровня моря и генерируется цветовая карта с объектом границы цвета. Например, для уровня моря 3000 м, где каждые 10 м отображаются цветом, нужно увеличить ocean на 300 и уменьшить terrain на 300 цветов:

def change_sea_level(new_sea_level):
    new_sea_level = new_sea_level / 10
    undersea = int(800 + new_sea_level)
    land = int(1400 - new_sea_level)
    colors_undersea = plt.cm.ocean(np.linspace(0.2, 0.8, undersea))
    undersea_map = LinearSegmentedColormap.from_list('undersea_map', colors_undersea, N=undersea)

    colors_land = plt.cm.terrain(np.linspace(0.25, 1, land))
    land_map = LinearSegmentedColormap.from_list('land_map', colors_land, N=land)

    colors = np.vstack((colors_undersea, colors_land))
    terrain_map = LinearSegmentedColormap.from_list('cut_terrain', colors, N=2221)

    bounds = np.arange(-8210, 14000, 10)
    norm = BoundaryNorm(bounds, terrain_map.N)
    return terrain_map, norm


terrain_map0, norm0 = change_sea_level(new_sea_level = 0)
terrain_map3000, norm3000 = change_sea_level(new_sea_level = 3000)
terrain_map_3000, norm_3000 = change_sea_level(new_sea_level = -3000)

Посмотрим, как выглядел бы Марс при уровне моря на 1500 м выше средней высоты. Такой Марс на 70% покрыт водой, как и Земля:

terrain_map1500, norm1500 = change_sea_level(new_sea_level = 1500)

fig, ax = plt.subplots()
fig.set_size_inches(14, 7)

ax.imshow(mars[0], cmap=terrain_map1500, norm=norm1500)
ax = label_features_terrain(ax)
ax.axis('off')
newax = fig.add_axes([0.82, 0.13, 0.08, 0.08], anchor='NE')
newax.imshow(logo)
newax.axis('off')
txt = ax.text(0.0, 0.02, "Martian Topography \n@PythonMaps",
              size=4,
              color='white',
              transform = ax.transAxes)
plt.show()

Отмывка рельефа

Добавим рельефу немного света, чтобы придать карте эффект 3D. Отмывка рельефа — это представление поверхности в 3D, обычно в оттенках серого. Более тёмные и светлые цвета — это тени и освещённые участки, визуально обнаруживаемые на модели рельефа. Отмывку рельефа часто применяют на картах, чтобы подчеркнуть трёхмерность и сделать данные визуально интересными. 

Чтобы сгенерировать данные с отмывкой рельефа, воспользуемся функцией earthpy. Есть два параметра, при настройке которых в зависимости от набора данных получаются очень разные результаты:

• Первый — это значение azimuth (азимут), которое варьируется от 0 до 360° и определяет источник света. 0° соответствует источнику света, указывающему строго на север.

• Второй — это altitude (высота), на которой источник света находится. Его значения варьируются от 1 до 90:

import earthpy as et
import earthpy.spatial as es

hillshade = es.hillshade(mars[0], azimuth=250, altitude=1)
fig, ax = plt.subplots(facecolor='#FCF6F5FF')
fig.set_size_inches(14, 7)
ax.imshow(mars[0], cmap=terrain_map, norm=norm)
ax.imshow(hillshade, cmap="Greys", alpha=0.2)
ax = label_features_terrain(ax)
ax.axis('off')
newax = fig.add_axes([0.8, 0.78, 0.08, 0.08], anchor='NE')
newax.imshow(logo)
newax.axis('off')
txt = ax.text(0.02, 0.03, "Martian Topography \n@PythonMaps",
              size=8,
              color='black',
              transform = ax.transAxes)
plt.show()

Итог — та же карта, но с объектами крупнее, особенно в северных океанах:

Заключение

Вот и всё. Мы создали красивую карту с чёткой визуализацией данных Марса. Впереди много статей о том, как делать геопространственные данные такими крутыми с помощью Python. 

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

Ссылки

Полный список использованных материалов.

Продолжить погружение в Data Science или Python вы сможете на наших курсах:

• Профессия Data Scientist

• Профессия Fullstack-разработчик на Python

Узнайте подробности здесь.