Documentation for PISM, the Parallel Ice Sheet Model

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
doc_misc [2014/01/13 14:07]
Ed Bueler
doc_misc [2014/01/20 13:49] (current)
Constantine Khroulev
Line 9: Line 9:
 Below is an example plotting csurf (the magnitude of the horizontal ice velocity at the ice surface) using a logarithmic color scale and marking the 100 <​sup>​m</​sup>/<​sub>​year</​sub>​ contour in black: Below is an example plotting csurf (the magnitude of the horizontal ice velocity at the ice surface) using a logarithmic color scale and marking the 100 <​sup>​m</​sup>/<​sub>​year</​sub>​ contour in black:
  
-[[https://raw.github.com/​pism/​pismdocs-template/master/​public_scripts/​basemap_ex1.py|basemap_ex1.py]]+++++ Click here to see the source of basemap_ex1.py | 
 +<code python basemap_ex1.py>​ 
 +#​!/​usr/​bin/​env python 
 +from mpl_toolkits.basemap import Basemap, NetCDFFile 
 +import numpy as np 
 +import matplotlib.pyplot as plt 
 +from matplotlib import colors 
 + 
 +nc = NetCDFFile("​g10km_0.nc",​ '​r'​) 
 + 
 +# we need to know longitudes and latitudes corresponding to out grid 
 +lon = nc.variables['​lon'​][:
 +lat = nc.variables['​lat'​][:​] 
 + 
 +# x and y *in the dataset* are only used to determine plotting domain 
 +# dimensions 
 +x = nc.variables['​x'​][:​] 
 +y = nc.variables['​y'​][:​] 
 +width = x.max() ​x.min() 
 +height = y.max() - y.min() 
 + 
 +# load data and mask out ice-free areas 
 +fill  = nc.variables['​csurf'​]._FillValue 
 +csurf = np.squeeze(nc.variables['​csurf'​][:​]) 
 +csurf = np.ma.array(csurf,​ mask = csurf == fill) 
 + 
 +m = Basemap(width=2*width, ​     # width in projection coordinates,​ in meters 
 +            height=height, ​     # height 
 +            resolution='​l', ​    # coastline resolution, can be '​l'​ (low), '​h'​ 
 +                                # (high) and '​f'​ (full) 
 +            projection='​stere',​ # stereographic projection 
 +            lat_ts=71, ​         # latitude of true scale 
 +            lon_0=-41, ​         # longitude of the plotting domain center 
 +            lat_0=72) ​          # latitude of the plotting domain center 
 + 
 +#​m.drawcoastlines() 
 + 
 +# draw the Blue Marble background (requires PIL, the Python Imaging Library) 
 +m.bluemarble() 
 + 
 +# convert longitudes and latitudes to x and y: 
 +xx,yy = m(lon, lat) 
 + 
 +# mark 100 m/a contours in black: 
 +m.contour(xx, yy, csurf, [100], colors="​black"​) 
 + 
 +# plot csurf using log color scale 
 +m.pcolormesh(xx,​yy,​csurf,​ 
 +             ​norm=colors.LogNorm(vmin=1,​ vmax=5e3)) # use log color scale, 
 +                                                     # omit this to use linear 
 +                                                     # color scale 
 + 
 +# add a colorbar: 
 +plt.colorbar(extend='​both',​ 
 +             ​ticks=[2,​ 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 3000]
 +             ​format="​%d"​) 
 + 
 +# draw parallels and meridians. The labels argument specifies where to draw 
 +# ticks: [left, right, top, bottom] 
 +m.drawparallels(np.arange(-55.,​90.,​5.),​ labels = [1, 0, 0, 0]) 
 +m.drawmeridians(np.arange(-120.,​30.,​10.),​ labels = [0, 0, 0, 1]) 
 + 
 +plt.title(r"​Modeled ice surface velocity, m/​year"​) 
 +plt.savefig('​g10km_0_csurf.png'​) 
 +</​code>​ 
 +++++ 
  
 The following alternative python code does nearly the same thing but uses [[https://​code.google.com/​p/​netcdf4-python/​|netcdf4-python]] to read NetCDF: The following alternative python code does nearly the same thing but uses [[https://​code.google.com/​p/​netcdf4-python/​|netcdf4-python]] to read NetCDF:
  
-[[https://raw.github.com/​pism/​pismdocs-template/master/​public_scripts/​basemap_ex2.py|basemap_ex2.py]]+++++ Click here to see the source of basemap_ex2.py | 
 +<code python basemap_ex2.py>​ 
 +#​!/​usr/​bin/​env python 
 +from mpl_toolkits.basemap import Basemap 
 + 
 +from netCDF4 import Dataset as NC 
 + 
 +import numpy as np 
 +import matplotlib.pyplot as plt 
 +from matplotlib import colors 
 +import sys 
 + 
 +rootname='​g10km_0'​ 
 + 
 +try: 
 +    nc = NC(rootname+'​.nc',​ '​r'​) 
 +except: 
 +    print "​ERROR:​ can't read from file ..." 
 +    sys.exit(1) 
 + 
 +# we need to know longitudes and latitudes corresponding to out grid 
 +lon = nc.variables['​lon'​][:
 +lat = nc.variables['​lat'​][:​] 
 + 
 +# x and y *in the dataset* are only used to determine plotting domain 
 +# dimensions 
 +x = nc.variables['​x'​][:​] 
 +y = nc.variables['​y'​][:​] 
 +width = x.max() ​x.min() 
 +height = y.max() - y.min() 
 + 
 +# load data and mask out ice-free areas 
 +fill  = nc.variables['​csurf'​]._FillValue 
 +csurf = np.squeeze(nc.variables['​csurf'​][:​]) 
 +csurf = np.ma.array(csurf,​ mask = csurf == fill) 
 + 
 +m = Basemap(width=1.2*width, ​     # width in projection coordinates,​ in meters 
 +            height=height, ​     # height 
 +            resolution='​l', ​    # coastline resolution, can be '​l'​ (low), '​h'​ 
 +                                # (high) and '​f'​ (full) 
 +            projection='​stere',​ # stereographic projection 
 +            lat_ts=71, ​         # latitude of true scale 
 +            lon_0=-41, ​         # longitude of the plotting domain center 
 +            lat_0=72) ​          # latitude of the plotting domain center 
 + 
 +#​m.drawcoastlines() 
 + 
 +# draw the Blue Marble background (requires PIL, the Python Imaging Library) 
 +m.bluemarble() 
 + 
 +# convert longitudes and latitudes to x and y: 
 +xx,yy = m(lon, lat) 
 + 
 +# mark 100 m/a contour in black: 
 +m.contour(xx, yy, csurf, [100], colors="​black"​) 
 + 
 +# plot csurf using log color scale 
 +m.pcolormesh(xx,​yy,​csurf,​ 
 +             ​norm=colors.LogNorm(vmin=1,​ vmax=6e3)) # use log color scale, 
 +                                                     # omit this to use linear 
 +                                                     # color scale 
 + 
 +# add a colorbar: 
 +plt.colorbar(extend='​both',​ 
 +             ​ticks=[2,​ 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000]
 +             ​format="​%d"​) 
 + 
 +# draw parallels and meridians. The labels argument specifies where to draw 
 +# ticks: [left, right, top, bottom] 
 +m.drawparallels(np.arange(-55.,​90.,​5.),​ labels = [1, 0, 0, 0]) 
 +m.drawmeridians(np.arange(-120.,​30.,​10.),​ labels = [0, 0, 0, 1]) 
 + 
 +plt.title(r"​modeled surface velocity, m/​year"​) 
 +plt.savefig(rootname+'​_csurf.png'​) 
 +</​code>​ 
 +++++ 
  
 This is the output of basemap_ex1.py:​ This is the output of basemap_ex1.py:​
 {{ g10km_0_csurf.png?​500 }} {{ g10km_0_csurf.png?​500 }}
doc_misc.txt · Last modified: 2014/01/20 13:49 by Constantine Khroulev
© 2016 by PISM | webmaster