I'm trying to create a map of Malawi with altitude shown. Something like this, but of Malawi of course:
I have downloaded some elevation data from here: http://research.jisao.washington.edu/data_sets/elevation/
This is a print of that data after I created a cube:
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
I started with importing my data, forming a cube, removing the extra variables (time and history) and limiting my data to the latitudes and longitudes for Malawi.
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')
This works well and the output is as follows:
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'})
However when I try to plot it, it doesn't work... this is what I did:
#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()
The error I get is 'Attribute Error: Unknown property type cmap'
and the following map of the whole world...
Any ideas?
I'll prepare the data the same as you, except to remove the
time
dimension I'll useiris.util.squeeze
, which removes any length-1 dimension.As @ImportanceOfBeingErnest says, you want a contour plot. When unsure what plotting function to use, I recommend browsing the matplotlib gallery to find something that looks similar to what you want to produce. Click on an image and it shows you the code.
So, to make the contour plot you can use the
matplotlib.pyplot.contourf
function, but you have to get the relevant data from the cube in the form ofnumpy
arrays:However,
iris
provides a shortcut to themaplotlib.pyplot
functions in the form ofiris.plot
. This automatically sets up an axes instance with the right projection, and passes the data from the cube through tomatplotlib.pyplot
. So the last two lines can simply become:There is also
iris.quickplot
, which is basically the same asiris.plot
, except that it automatically adds a colorbar and labels where appropriate:Once plotted, you can get hold of the axes instance and add your other items (for which I simply copied your code):