python-2.7matplotlibcartopypython-iris

Plotting Elevation in Python


I'm trying to create a map of Malawi with altitude shown. Something like this, but of Malawi of course:

enter image description here

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...

enter image description here

Any ideas?


Solution

  • 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)
    

    matplotlib.pyplot version

    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)
    

    iris.plot version

    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)
    

    iris.quickplot version

    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)
    

    iris.quickplot with additions