在Python中绘制高程

Question

我正在尝试创建一张显示高度的马拉维地图。这样的事情，但马拉维当然​​：

我从这里下载了一些高程数据：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'和以下全世界地图......

有任何想法吗？

Answer 1

我将准备与您相同的数据，除了删除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 gallery以找到与您想要生成的内容类似的内容。单击图像，它会显示代码。

因此，要制作轮廓图，您可以使用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)

然而，iris以maplotlib.pyplot的形式提供了iris.plot函数的捷径。这会使用正确的投影自动设置轴实例，并将数据从立方体传递到matplotlib.pyplot。所以最后两行可以简单地变成：

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

还有iris.quickplot，它与iris.plot基本相同，只是它会在适当时自动添加颜色条和标签：

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

绘制完成后，您可以获取轴实例并添加其他项目（我只是复制了您的代码）：

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)

Answer 2

看起来你想要一个类似轮廓图的东西。而不是

plot = ax.plot(...)

你可能想用

plot = ax.contourf(...)

很可能你也想把经度和经度作为contourf的参数，

plot = ax.contourf(longitude, latitude, Elev, ...)

Answer 3

你可以尝试添加这个：

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