This is part of Python for Geosciences notes.
=============
Some of the things will not work on ZMAW computers (Iris, Cartopy).
Iris seeks to provide a powerful, easy to use, and community-driven Python library for analysing and visualising meteorological and oceanographic data sets. Kind of Ferret replacement. Developed in the Met Office by group of 7 full time developers. There are more than 300 active python users in Met Office.
With Iris you can:
Here we load data from netCDF file in to cube object and plot first time step. Note automatic unpacking, and the title.
import iris
import iris.quickplot as qplt
temperature = iris.load_cube('air.sig995.2012.nc')
qplt.contourf(temperature[0,:,:])
gca().coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x4a05950>
This is how iris cube look like:
print temperature
4xDaily Air temperature at sigma level 995 / (degK) (time: 1464; latitude: 73; longitude: 144) Dimension coordinates: time x - - latitude - x - longitude - - x Attributes: Conventions: COARDS GRIB_id: 11 GRIB_name: TMP actual_range: [ 191.1000061 323. ] dataset: NMC Reanalysis description: Data is from NMC initialized reanalysis (4x/day). These are the 0.9950... history: created 2011/12 by Hoop (netCDF2.3) least_significant_digit: 1 level_desc: Surface parent_stat: Other platform: Model precision: 2 references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html statistic: Individual Obs title: 4x daily NMC reanalysis (2012) unpacked_valid_range: [ 185.16000366 331.16000366] var_desc: Air temperature
We can perform different operations on cubes. For example create zonal mean:
zonal_mean = temperature.collapsed('latitude', iris.analysis.MEAN)
/usr/local/lib/python2.7/dist-packages/Iris-1.6.0_dev-py2.7.egg/iris/coords.py:855: UserWarning: Collapsing a non-contiguous coordinate. Metadata may not be fully descriptive for "latitude". self.name())
qplt.contourf(zonal_mean)
<matplotlib.contour.QuadContourSet instance at 0x4a14560>
Here we plot timesiries from one point:
#Code is a bit more complicated in order to fix issue with dates formating
fig = figure()
qplt.plot(temperature[:,10,10])
fig.autofmt_xdate()
Section along longitude:
qplt.plot(temperature[0,:,10])
[<matplotlib.lines.Line2D at 0x7255b90>]
Cartopy is a library providing cartographic tools for python.
Some of the key features of cartopy are:
Simple plot:
import cartopy.crs as ccrs
ax = plt.axes(projection=ccrs.PlateCarree())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x75d8290>
We only change projection:
ax = plt.axes(projection=ccrs.Mollweide())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x77d6550>
ax = plt.axes(projection=ccrs.Robinson())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7c3b190>
fig = figure(figsize=(7,7))
ax = plt.axes(projection=ccrs.NorthPolarStereo())
ax.set_extent([0, 360, 50, 90], crs=ccrs.PlateCarree())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x80d6c90>
One of the things it was originally created for is to handle datelines properly. Below is an example form the post by Filipe Fernandes, that demonstrate this feature. He show how to plot Challenger Expedition track in Basemap and Cartopy.
First download the data:
!wget https://raw.github.com/ocefpaf/ocefpaf.github.io/master/downloads/notebooks/data/challenger_path.csv
kw = dict(color='#FF9900', linestyle='-', linewidth=1.5)
lon, lat = np.loadtxt('./challenger_path.csv', delimiter=',', unpack=True)
from mpl_toolkits.basemap import Basemap
def make_basemap(projection='robin', figsize=(10, 5), resolution='c'):
fig, ax = plt.subplots(figsize=figsize)
m = Basemap(projection=projection, resolution=resolution,
lon_0=0, ax=ax)
m.drawcoastlines()
m.fillcontinents(color='0.85')
parallels = np.arange(-60, 90, 30.)
meridians = np.arange(-360, 360, 60.)
m.drawparallels(parallels, labels=[1, 0, 0, 0])
m.drawmeridians(meridians, labels=[0, 0, 1, 0])
return fig, m
fig, m = make_basemap()
_ = m.plot(*m(lon, lat), **kw)
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def make_cartopy(projection=ccrs.Robinson(), figsize=(10, 5), resolution='110m'):
fig, ax = plt.subplots(figsize=figsize, subplot_kw=dict(projection=projection))
ax.set_global()
ax.coastlines(resolution=resolution, color='k')
gl = ax.gridlines(draw_labels=False) # Only PlateCarree and Mercator plots are currently supported.
ax.add_feature(cfeature.LAND, facecolor='0.75')
return fig, ax
fig, ax = make_cartopy(projection=ccrs.Robinson(), resolution='110m')
_ = ax.plot(lon, lat, transform=ccrs.Geodetic(), **kw)
Pydap is a pure Python library implementing the Data Access Protocol, also known as DODS or OPeNDAP. You can use Pydap as a client to access hundreds of scientific datasets in a transparent and efficient way through the internet; or as a server to easily distribute your data from a variety of formats.
from pydap.client import open_url
We are going to access sea ice data from CliSAP-Integrated Climate Data Center (ICDC)
dataset = open_url("http://icdc.zmaw.de/thredds/dodsC/amsre_asi_nh_2011")
print dataset
{'latitude': <pydap.model.BaseType object at 0x932dd90>, 'longitude': <pydap.model.BaseType object at 0x932ddd0>, 'time': <pydap.model.BaseType object at 0x91dea50>, 'icecon': <pydap.model.BaseType object at 0x91dead0>}
ice = dataset['icecon']
ice.shape
(273, 1792, 1216)
ice.attributes
{'_FillValue': -1000, 'missing_value': -1000, 'scale_factor': 0.10000000149011612, 'units': '%'}
imshow(squeeze(ice[0,:,:])*0.10000000149011612)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x8670368>
The F2PY project is created to unify the efforts of supporting easy connection between Fortran and Python languages. Example below is from Using Python and FORTRAN with F2py.
Create FORTRAN file (use %%file instead of %%writefile if you on IPython < 1.0):
%%writefile FIB1.F
C FILE: FIB1.F
SUBROUTINE FIB(A,N)
C
C CALCULATE FIRST N FIBONACCI NUMBERS
C
INTEGER N
REAL*8 A(N)
DO I=1,N
IF (I.EQ.1) THEN
A(I) = 0.0D0
ELSEIF (I.EQ.2) THEN
A(I) = 1.0D0
ELSE
A(I) = A(I-1) + A(I-2)
ENDIF
ENDDO
END
C END FILE FIB1.F
Overwriting FIB1.F
Compile it with f2py:
!f2py -c -m --fcompiler=gnu95 fib1 FIB1.F
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "fib1" sources
f2py options: []
f2py:> /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.c
creating /tmp/tmpId9Fu8
creating /tmp/tmpId9Fu8/src.linux-x86_64-2.7
Reading fortran codes...
Reading file 'FIB1.F' (format:fix,strict)
Post-processing...
Block: fib1
Block: fib
Post-processing (stage 2)...
Building modules...
Building module "fib1"...
Constructing wrapper function "fib"...
fib(a,[n])
Wrote C/API module "fib1" to file "/tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.c"
adding '/tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.c' to sources.
adding '/tmp/tmpId9Fu8/src.linux-x86_64-2.7' to include_dirs.
copying /usr/local/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpId9Fu8/src.linux-x86_64-2.7
copying /usr/local/lib/python2.7/dist-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpId9Fu8/src.linux-x86_64-2.7
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
customize Gnu95FCompiler
Found executable /usr/bin/gfortran
customize Gnu95FCompiler using build_ext
building 'fib1' extension
compiling C sources
C compiler: gcc -pthread -fno-strict-aliasing -DNDEBUG -g -fwrapv -O2 -Wall -Wstrict-prototypes -fPIC
creating /tmp/tmpId9Fu8/tmp
creating /tmp/tmpId9Fu8/tmp/tmpId9Fu8
creating /tmp/tmpId9Fu8/tmp/tmpId9Fu8/src.linux-x86_64-2.7
compile options: '-I/tmp/tmpId9Fu8/src.linux-x86_64-2.7 -I/usr/local/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'
gcc: /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.c
In file included from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1728:0,
from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,
from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:15,
from /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.h:13,
from /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.c:18:
/usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
/tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.c:111:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]
gcc: /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.c
In file included from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarraytypes.h:1728:0,
from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/ndarrayobject.h:17,
from /usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/arrayobject.h:15,
from /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.h:13,
from /tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.c:2:
/usr/local/lib/python2.7/dist-packages/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
compiling Fortran sources
Fortran f77 compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /usr/bin/gfortran -Wall -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -Wall -fno-second-underscore -fPIC -O3 -funroll-loops
compile options: '-I/tmp/tmpId9Fu8/src.linux-x86_64-2.7 -I/usr/local/lib/python2.7/dist-packages/numpy/core/include -I/usr/include/python2.7 -c'
gfortran:f77: FIB1.F
/usr/bin/gfortran -Wall -Wall -shared /tmp/tmpId9Fu8/tmp/tmpId9Fu8/src.linux-x86_64-2.7/fib1module.o /tmp/tmpId9Fu8/tmp/tmpId9Fu8/src.linux-x86_64-2.7/fortranobject.o /tmp/tmpId9Fu8/FIB1.o -lgfortran -o ./fib1.so
running scons
Removing build directory /tmp/tmpId9Fu8
Import resulting fib1.so as python library:
import fib1
Read some auto generated documentation:
print fib1.__doc__
print fib1.fib.__doc__
This module 'fib1' is auto-generated with f2py (version:2). Functions: fib(a,n=len(a)) . fib - Function signature: fib(a,[n]) Required arguments: a : input rank-1 array('d') with bounds (n) Optional arguments: n := len(a) input int
Use fib function:
a=zeros(15,'d')
fib1.fib(a)
print a
[ 0. 1. 1. 2. 3. 5. 8. 13. 21. 34. 55. 89. 144. 233. 377.]
netCDF4-python is advanced Python interface to the netCDF version 4 library.
Among other features it can read data from a multi-file netCDF dataset.
Let's download one more file from NCEP reanalysis data:
!wget ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/surface/air.sig995.2011.nc
from netCDF4 import MFDataset
f = MFDataset('air.sig995.????.nc')
f.variables
OrderedDict([(u'lat', <netCDF4.Variable object at 0x91e22d0>), (u'lon', <netCDF4.Variable object at 0x91e2450>), (u'time', <netCDF4._Variable object at 0xa4e0790>), (u'air', <netCDF4._Variable object at 0xa4e0390>)])
air = f.variables['air']
time = f.variables['time']
lat = f.variables['lat']
The air variable now have 2 years of data:
air.shape
(2924, 73, 144)
It also have very nice functions for dates processing:
from netCDF4 import num2date
time.units
u'hours since 1-1-1 00:00:0.0'
We can convert our time values to the datetime format, that python can work with:
time_conv = num2date(time, time.units)
time_conv
array([datetime.datetime(2011, 1, 1, 0, 0), datetime.datetime(2011, 1, 1, 6, 0), datetime.datetime(2011, 1, 1, 12, 0), ..., datetime.datetime(2012, 12, 31, 6, 0), datetime.datetime(2012, 12, 31, 12, 0), datetime.datetime(2012, 12, 31, 18, 0)], dtype=object)
fig = figure()
plot(time_conv, air[:,10,10])
fig.autofmt_xdate()
Note that we don't have to apply scale factor and add offset, it's done automatically.
contourf(lat[:],time_conv,air[:,:,10])
colorbar()
<matplotlib.colorbar.Colorbar instance at 0xc90db90>