| 1 | import numpy as np |
|---|
| 2 | from netCDF4 import Dataset as NetCDFFile |
|---|
| 3 | import matplotlib as mpl |
|---|
| 4 | mpl.use('Agg') |
|---|
| 5 | import matplotlib.pyplot as plt |
|---|
| 6 | from mpl_toolkits.basemap import Basemap |
|---|
| 7 | |
|---|
| 8 | filename = '/media/ExtDiskC_ext4/bkup/foudre/feinetes/BootCamp/files/Relief.nc' |
|---|
| 9 | toponame = 'RELIEF' |
|---|
| 10 | ncobj = NetCDFFile(filename, 'r') |
|---|
| 11 | |
|---|
| 12 | objtopo = ncobj.variables[toponame] |
|---|
| 13 | olon = ncobj.variables['longitude'] |
|---|
| 14 | olat = ncobj.variables['latitude'] |
|---|
| 15 | topovals = objtopo[:] |
|---|
| 16 | utopo = objtopo.getncattr('units') |
|---|
| 17 | |
|---|
| 18 | plt.rc('text', usetex=True) |
|---|
| 19 | |
|---|
| 20 | dx = olon.shape[0] |
|---|
| 21 | dy = olat.shape[0] |
|---|
| 22 | nlon = np.min(olon) |
|---|
| 23 | xlon = np.max(olon) |
|---|
| 24 | nlat = np.min(olat) |
|---|
| 25 | xlat = np.max(olat) |
|---|
| 26 | plt.xlim(nlon, xlon) |
|---|
| 27 | plt.ylim(nlat, xlat) |
|---|
| 28 | m = Basemap(projection='cyl', \ |
|---|
| 29 | llcrnrlon=nlon, llcrnrlat=nlat, \ |
|---|
| 30 | urcrnrlon=xlon, urcrnrlat= xlat, \ |
|---|
| 31 | resolution='l') |
|---|
| 32 | |
|---|
| 33 | lons = olon[:] |
|---|
| 34 | lats = olat[:] |
|---|
| 35 | x,y = m(lons,lats) |
|---|
| 36 | |
|---|
| 37 | plt.pcolormesh(x,y,topovals, \ |
|---|
| 38 | cmap=plt.get_cmap('terrain'), \ |
|---|
| 39 | vmin=0., vmax=7000.) |
|---|
| 40 | cbar = plt.colorbar() |
|---|
| 41 | m.drawcoastlines() |
|---|
| 42 | cbar.set_label(utopo) |
|---|
| 43 | |
|---|
| 44 | plt.xlabel('W-E') |
|---|
| 45 | plt.ylabel('S-N') |
|---|
| 46 | plt.title('Topography from a ' + \ |
|---|
| 47 | 'LMDZ domain at 1/' + \ |
|---|
| 48 | str(olon.shape[0]/360) + ' deg') |
|---|
| 49 | |
|---|
| 50 | plt.savefig("topo.png") |
|---|
| 51 | plt.close() |
|---|
