Построение высоты в 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' и следующую карту всего мира...

введите здесь описание изображения

Любые идеи?


person ErikaAWT    schedule 19.12.2017    source источник


Ответы (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 с дополнениями

person RuthC    schedule 21.12.2017
comment
Спасибо @RuthC, это очень четкий и полезный ответ! С Новым Годом! - person ErikaAWT; 03.01.2018

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

plot = ax.plot(...)

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

plot = ax.contourf(...)

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

plot = ax.contourf(longitude, latitude, Elev, ...)
person ImportanceOfBeingErnest    schedule 19.12.2017
comment
Спасибо @ImportanceOfBeingErnest. Я изменил его на plot=ax.contourf(Elev.coord('latitude'), Elev.coord('longitude'), Elev.data, cmap=YlGn, levels=np.arange(0,2000,150), extend='both'), но теперь получаю сообщение об ошибке: ValueError: setting an array element with a sequence. - person ErikaAWT; 19.12.2017
comment
Всегда предоставляйте полный код и полную трассировку ошибок. Иначе невозможно узнать, что происходит. - person ImportanceOfBeingErnest; 19.12.2017
comment
Извините, я обновил #plot altitude data plot=ax.plot(Elev, cmap=mpl_cm.get_cmap('YlGn'), levels=np.arange(0,2000,150), extend='both'), чтобы прочитать это #plot altitude data YlGn = plt.get_cmap('YlGn') plot=ax.contourf(Elev.coord('latitude'), Elev.coord('longitude'), Elev.data, cmap=YlGn, levels=np.arange(0,2000,150), extend='both') - person ErikaAWT; 19.12.2017
comment
Ошибка: файл /home/s0899345/datastore/Climate_Modelling/Python_Code_and_Output_Images/Map_boundary_altitudes.py, строка 51, в main plot=ax.contourf(Elev.coord('широта'), Elev.coord('долгота'), Elev.coord('долгота'), Elev. data, cmap=YlGn, level=np.arange(0,2000,150), extend='both') Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py , строка 1333, в контуре результат = matplotlib.axes.Axes.contourf(self, *args, **kwargs) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py, строка 1819, функция внутреннего возврата (ax, *args, **kwargs) - person ErikaAWT; 19.12.2017
comment
(Продолжение) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/axes/_axes.py, строка 5627, в convertf возвращает mcontour.QuadContourSet(self, *args, **kwargs) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/contour.py, строка 1424, в init ContourSet.__init__(self, ax, *args, **kwargs ) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/contour.py, строка 863, в init self._process_args(*args, **kwargs) - person ErikaAWT; 19.12.2017
comment
(Продолжение 2) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/contour.py, строка 1445, в _process_args x, y, z = self._contour_args(args, kwargs) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/matplotlib/contour.py, строка 1528, в _contour_args x, y, z = self._check_xyz(args[:3], kwargs) Файл /scratch/ s0899345/anaconda/lib/python2.7/site-packages/matplotlib/contour.py, строка 1557, в _check_xyz x = np.asarray(x, dtype=np.float64) - person ErikaAWT; 19.12.2017
comment
(Продолжение 3) Файл /scratch/s0899345/anaconda/lib/python2.7/site-packages/numpy/core/numeric.py, строка 531, в массиве возврата asarray (a, dtype, copy=False, order= порядок) ValueError: установка элемента массива с последовательностью. - person ErikaAWT; 19.12.2017
comment
Бесполезно помещать код или трассировку в комментарии. В любом случае кажется, что вы не можете поместить объект DimCoords в контур, но вместо этого вам нужно напрямую передать массив numpy. Я никогда не работал с Айрис, но я бы попробовал Elev.coord('longitude').data. - person ImportanceOfBeingErnest; 19.12.2017
comment
Извините, @ImportanceOfBeingErnest. Еще раз спасибо, но это тоже не работает. Вы получаете ошибку AttributeError: 'DimCoord' object has no attribute 'data' Как вы можете видеть в моем исходном сообщении, я напечатал долготу и широту как Elev.coord('longitude'), и это работает. - person ErikaAWT; 19.12.2017
comment
Я попытался предопределить их как Lat = Elev.coord('latitude'), Lon = Elev.coord('longitude') и Elev = Elev.data, а затем запустить plot=ax.contourf(Lat, Lon, Elev, cmap=YlGn, levels=np.arange(0,2000,150), extend='both'), но это все равно возвращается с исходной ошибкой ValueError: setting an array element with a sequence. - person ErikaAWT; 19.12.2017
comment
Вам нужно выяснить, как получить массив numpy из объекта DimCoord. Я не могу помочь вам с этим. - person ImportanceOfBeingErnest; 19.12.2017

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

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 
person Vasily Bronsky    schedule 19.12.2017
comment
Спасибо Василий, Спасибо за вашу помощь, но это не сработало :( Без изменений. - person ErikaAWT; 19.12.2017
comment
вы можете попытаться найти что-нибудь здесь - person Vasily Bronsky; 19.12.2017
comment
Спасибо Василий, я там уже смотрел, но безрезультатно. Это странно, поскольку я использовал ту же кодировку в другом наборе данных, и она отлично работает. Мне не нужно было ничего конкретно импортировать, и cmap понятен. Я просто не могу заставить его работать с этими данными.... - person ErikaAWT; 19.12.2017