在 Python 中绘制高程
Plotting Elevation in 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'和下面的全世界地图...
有什么想法吗？
你可以试试加上这个：
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
您似乎想要等高线图之类的东西。所以而不是
plot = ax.plot(...)
您可能想使用
plot = ax.contourf(...)
您很可能还想将纬度和经度作为参数提供给 contourf、
plot = ax.contourf(longitude, latitude, Elev, ...)
除了要删除 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 以 iris.plot 的形式提供了 maplotlib.pyplot 函数的快捷方式。这会自动设置一个具有正确投影的轴实例，并将数据从立方体传递到 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)
我正在尝试创建一张显示海拔高度的马拉维地图。像这样，但当然是马拉维：
我从这里下载了一些海拔数据：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'和下面的全世界地图...
有什么想法吗？
你可以试试加上这个：
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
您似乎想要等高线图之类的东西。所以而不是
plot = ax.plot(...)
您可能想使用
plot = ax.contourf(...)
您很可能还想将纬度和经度作为参数提供给 contourf、
plot = ax.contourf(longitude, latitude, Elev, ...)
除了要删除 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 以 iris.plot 的形式提供了 maplotlib.pyplot 函数的快捷方式。这会自动设置一个具有正确投影的轴实例，并将数据从立方体传递到 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)