Возведение высот в Python

Я пытаюсь создать карту Малави с показанной высотой. Как то так, но из Малави конечно

Я скачал некоторые данные по высоте здесь: http://research.jisao.washington.edu/data_sets/elevation/

Это распечатка этих данных после того, как я создал куб:

meters, from 5-min data / (unknown) (time: 1; latitude: 360; longitude: 720)
     Dimension coordinates:
          time                           x            -               -
          latitude                       -            x               -
          longitude                      -            -               x
     Attributes:
          history: 
Elevations calculated from the TBASE 5-minute
latitude-longitude resolution...
          invalid_units: meters, from 5-min data

Я начал с импорта моих данных, формирования куба, удаления дополнительных переменных (времени и истории) и ограничения моих данных широтами и долготами Малави.

import matplotlib.pyplot as plt
import matplotlib.cm as mpl_cm
import numpy as np
import iris
import cartopy
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import iris.analysis.cartography


def main():

    #bring in altitude data
    Elev = '/exports/csce/datastore/geos/users/s0899345/Climate_Modelling/Actual_Data/elev.0.5-deg.nc'

    Elev= iris.load_cube(Elev)

    #remove variable for time
    del Elev.attributes['history']
    Elev = Elev.collapsed('time', iris.analysis.MEAN)

    Malawi = iris.Constraint(longitude=lambda v: 32.0 <= v <= 36., latitude=lambda v: -17. <= v <= -8.)      
    Elev = Elev.extract(Malawi)

    print 'Elevation'
    print Elev.data
    print 'latitude'
    print Elev.coord('latitude')
    print 'longitude'
    print Elev.coord('longitude')

Это работает хорошо, и результат выглядит следующим образом:

Elevation
[[  978.  1000.  1408.  1324.  1080.  1370.  1857.  1584.]
 [ 1297.  1193.  1452.  1611.  1354.  1480.  1350.   627.]
 [ 1418.  1490.  1625.  1486.  1977.  1802.  1226.   482.]
 [ 1336.  1326.  1405.   728.  1105.  1559.  1139.   789.]
 [ 1368.  1301.  1463.  1389.   671.   942.   947.   970.]
 [ 1279.  1116.  1323.  1587.   839.  1014.  1071.  1003.]
 [ 1096.   969.  1179.  1246.   855.   979.   927.   638.]
 [  911.   982.  1235.  1324.   681.   813.   814.   707.]
 [  749.   957.  1220.  1198.   613.   688.   832.   858.]
 [  707.  1049.  1037.   907.   624.   771.  1142.  1104.]
 [  836.  1044.  1124.  1120.   682.   711.  1126.   922.]
 [ 1050.  1204.  1199.  1161.   777.   569.   999.   828.]
 [ 1006.   869.  1183.  1230.  1354.   616.   762.   784.]
 [  838.   607.   883.  1181.  1174.   927.   591.   856.]
 [  561.   402.   626.   775.  1053.   726.   828.   733.]
 [  370.   388.   363.   422.   508.   471.   906.  1104.]
 [  504.   326.   298.   208.   246.   160.   458.   682.]
 [  658.   512.   334.   309.   156.   162.   123.   340.]]
latitude
DimCoord(array([ -8.25,  -8.75,  -9.25,  -9.75, -10.25, -10.75, -11.25, -11.75,
       -12.25, -12.75, -13.25, -13.75, -14.25, -14.75, -15.25, -15.75,
       -16.25, -16.75], dtype=float32), standard_name='latitude', units=Unit('degrees'), var_name='lat', attributes={'title': 'Latitude'})
longitude
DimCoord(array([ 32.25,  32.75,  33.25,  33.75,  34.25,  34.75,  35.25,  35.75], dtype=float32), standard_name='longitude', units=Unit('degrees'), var_name='lon', attributes={'title': 'Longitude'})

Однако, когда я пытаюсь построить это, это не работает... вот что я сделал:

#plot map with physical features 
ax = plt.axes(projection=cartopy.crs.PlateCarree())
ax.add_feature(cartopy.feature.COASTLINE)   
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
ax.add_feature(cartopy.feature.RIVERS)
#plot altitude data
plot=ax.plot(Elev, cmap=mpl_cm.get_cmap('YlGn'), levels=np.arange(0,2000,150), extend='both') 
#add colour bar index and a label
plt.colorbar(plot, label='meters above sea level')
#set map boundary
ax.set_extent([32., 36., -8, -17]) 
#set axis tick marks
ax.set_xticks([33, 34, 35]) 
ax.set_yticks([-10, -12, -14, -16]) 
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#save the image of the graph and include full legend
plt.savefig('Map_data_boundary', bbox_inches='tight')
plt.show()

Я получаю ошибку 'Attribute Error: Unknown property type cmap' и следующая карта всего мира...

Есть идеи?

Источник

3 ответа

Решение

Я подготовлю данные так же, как и вы, за исключением удаления time измерение я буду использовать iris.util.squeeze, который удаляет любое измерение длины-1.

import iris

elev = iris.load_cube('elev.0.5-deg.nc')
elev = iris.util.squeeze(elev)
malawi = iris.Constraint(longitude=lambda v: 32.0 <= v <= 36.,
                         latitude=lambda v: -17. <= v <= -8.)      
elev = elev.extract(malawi)

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

Итак, для построения контурного графика вы можете использовать matplotlib.pyplot.contourf функция, но вы должны получить соответствующие данные из куба в виде numpy массивы:

import matplotlib.pyplot as plt
import matplotlib.cm as mpl_cm
import numpy as np
import cartopy

cmap = mpl_cm.get_cmap('YlGn')
levels = np.arange(0,2000,150)
extend = 'max'

ax = plt.axes(projection=cartopy.crs.PlateCarree())
plt.contourf(elev.coord('longitude').points, elev.coord('latitude').points, 
             elev.data, cmap=cmap, levels=levels, extend=extend)

версия matplotlib.pyplot

Тем не мение, iris обеспечивает ярлык для maplotlib.pyplot функции в виде iris.plot, Это автоматически устанавливает экземпляр оси с правильной проекцией и передает данные из куба в matplotlib.pyplot, Таким образом, последние две строки могут просто стать:

import iris.plot as iplt
iplt.contourf(elev, cmap=cmap, levels=levels, extend=extend)

версия iris.plot

Существует также iris.quickplot что в основном так же, как iris.plot за исключением того, что он автоматически добавляет цветовую панель и метки, где это необходимо:

import iris.quickplot as qplt
qplt.contourf(elev, cmap=cmap, levels=levels, extend=extend)

версия iris.quickplot

После построения вы можете получить экземпляр оси и добавить другие элементы (для которых я просто скопировал ваш код):

from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

qplt.contourf(elev, cmap=cmap, levels=levels, extend=extend)
ax = plt.gca()

ax.add_feature(cartopy.feature.COASTLINE)   
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
ax.add_feature(cartopy.feature.RIVERS)

ax.set_xticks([33, 34, 35]) 
ax.set_yticks([-10, -12, -14, -16]) 
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

iris.quickplot с дополнениями

Источник

Кажется, вы хотите что-то вроде контурного сюжета. Так что вместо

plot = ax.plot(...)

вы, вероятно, хотите использовать

plot = ax.contourf(...)

Скорее всего, вы также хотите указать широту и долготу в качестве аргументов contourf,

plot = ax.contourf(longitude, latitude, Elev, ...)
Источник

Вы можете попробовать добавить это:

import matplotlib.colors as colors

color = plt.get_cmap('YlGn')  # and change cmap=mpl_cm.get_cmap('YlGn') to  cmap=color

А также попробуйте обновить ваш matplotlib:

pip install --upgrade matplotlib

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

color = plt.get_cmap('YlGn')  # and change cmap=mpl_cm.get_cmap('YlGn') to  cmap=color
Источник
Другие вопросы по тегам