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 use iris.util.squeeze
, which removes any length-1 dimension.
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)
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 of numpy
arrays:
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)
However, iris
provides a shortcut to the maplotlib.pyplot
functions in the form of iris.plot
. This automatically sets up an axes instance with the right projection, and passes the data from the cube through to matplotlib.pyplot
. So the last two lines can simply become:
import iris.plot as iplt
iplt.contourf(elev, cmap=cmap, levels=levels, extend=extend)
There is also iris.quickplot
, which is basically the same as iris.plot
, except that it automatically adds a colorbar and labels where appropriate:
import iris.quickplot as qplt
qplt.contourf(elev, cmap=cmap, levels=levels, extend=extend)
Once plotted, you can get hold of the axes instance and add your other items (for which I simply copied your code):
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)