代码之家 › 专栏 › 技术社区 › ErikaAWT

在Python中打印高程

python-iris cartopy matplotlib python-2.7

ErikaAWT · 技术社区 · 8 年前

我正试图绘制一张显示了海拔高度的马拉维地图。像这样的,当然是马拉维:

我从这里下载了一些高程数据: 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 回复 | 直到 8 年前

RuthC 8 年前

我将与您一样准备数据,只是要删除 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 ,基本上与 艾里斯。情节 ,但它会在适当的情况下自动添加颜色栏和标签:

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

绘制完成后,您可以获得axes实例并添加其他项目(我只是复制了您的代码):

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)

ImportanceOfBeingErnest 8 年前

看起来你想要一个等高线图。因此

plot = ax.plot(...)

您可能想使用

plot = ax.contourf(...)

很可能您还希望将纬度和经度作为 contourf ,

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

Vasily Bronsky 8 年前

您可以尝试添加以下内容:

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