This is part of Python for Geosciences notes.
=============
Some of the things will not work on ZMAW computers (Iris, Cartopy).
Iris¶
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:
- Use a single API to work on your data, irrespective of its original format.
- Read and write (CF-)netCDF, GRIB, and PP files.
- Easily produce graphs and maps via integration with matplotlib and cartopy.
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()
This is how iris cube look like:
print temperature
We can perform different operations on cubes. For example create zonal mean:
zonal_mean = temperature.collapsed('latitude', iris.analysis.MEAN)
qplt.contourf(zonal_mean)
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])
Cartopy¶
Cartopy is a library providing cartographic tools for python.
Some of the key features of cartopy are:
- object oriented projection definitions
- point, line, polygon and image transformations between projections
- integration to expose advanced mapping in matplotlib with a simple and intuitive interface
Simple plot:
import cartopy.crs as ccrs
ax = plt.axes(projection=ccrs.PlateCarree())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
We only change projection:
ax = plt.axes(projection=ccrs.Mollweide())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
ax = plt.axes(projection=ccrs.Robinson())
qplt.contourf(temperature[0,:,:])
gca().coastlines()
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()
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
Basemap version:¶
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)
Cartopy version:¶
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¶
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
ice = dataset['icecon']
ice.shape
ice.attributes
imshow(squeeze(ice[0,:,:])*0.10000000149011612)
colorbar()
F2PY¶
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
Compile it with f2py:
!f2py -c -m --fcompiler=gnu95 fib1 FIB1.F
Import resulting fib1.so as python library:
import fib1
Read some auto generated documentation:
print fib1.__doc__
print fib1.fib.__doc__
Use fib function:
a=zeros(15,'d')
fib1.fib(a)
print a
netCDF4-python¶
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
air = f.variables['air']
time = f.variables['time']
lat = f.variables['lat']
The air variable now have 2 years of data:
air.shape
It also have very nice functions for dates processing:
from netCDF4 import num2date
time.units
We can convert our time values to the datetime format, that python can work with:
time_conv = num2date(time, time.units)
time_conv
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()
Comments !