<<TableOfContents(4)>>

= PyNIO =


Find full PyNio documentation [[ http://www.pyngl.ucar.edu/Nio.shtml | here ]]

 * Download a [[ http://www.unidata.ucar.edu/software/netcdf/examples/sresa1b_ncar_ccsm3_0_run1_200001.nc | sample file (CCM precipitation flux, air temperature, etc)]]
 * load python module and start ipython 

{{{
module load python/2.7-ve0
ipython
}}}

Ipython listing


== Read NetCDF files ==

{{{#!python
# load Nio module
import Nio

# open netcdf file
f = Nio.open_file('sresa1b_ncar_ccsm3_0_run1_200001.nc', 'r') # 'r' stands for "read rights"

# check contents
print f

# get dimensions
dimNames = f.dimensions.keys()

# get the size of dimension
dimSize = f.dimensions['dimName']

# read single variable:
pr = f.variables['pr']
print pr

# read variables contents
pr_data = pr[:]
print pr_data

# print slice of the variable data
slice = pr_data[0,::-1,:]

# close file
f.close()
}}}

== Create NetCDF file ==

{{{#!python
# create new file
import Nio
import numpy

f = Nio.open_file('test.nc', 'w') # "w" stands for writing rights

# create a dimension
f.create_dimension('time',100)

# create variable
f.create_variable('temperature', 'i', ('time',)) # dimension 'time' must be created prior to this step

# dimensions and variables can have same names
f.create_variable('time', 'i', ('time',))

# Dimension types:
#    'd': 64 bit float
#    'f': 32 bit float
#    'l': long
#    'i': 32 bit integer
#    'h': 16 bit integer
#    'b': 8 bit integer
#    'S1': character 

# Create an attribute 
f.variables['temperature'].units = "K"

# Add scaling factor and an offset
# when the variable will be read, it will need to be multiplied by the scale_factor first and the added to the offset value
f.variables['temperature'].scale_factor = 0.1
f.variables['temperature'].add_offset = 273

# create variable contents and assign it to the NetCDF variable
temp = numpy.arange(0,100,1, dtype = numpy.int32)
f.variables['temperature'][:] = temp

time = numpy.arange(1000,1100,1, dtype = numpy.int32)
f.variables['time'][:] = time

f.close()

}}}

Any netcdf file can be reopened later and new dimensions and variables can be added.

== Install PyNio locally ==

Not in repositories, needs to be installed manually
web: http://www.pyngl.ucar.edu/Nio.shtml

requires: python-numpy (1.4.1, 1.5.1 or 1.6.0)

how to get:
• go to http://earthsystemgrid.org
• login with your open id (or request one)
• go to "software" directory and download PyNio 1.4.1 precompiled binaries
• pay attention which binaries will fit your system: Linux Debian or Redhat; MacOsX 10.6; there are different versions of gcc and numpy

install: tar -xzvf "package.tar.gz" and sudo python setup.py install

On MPI servers Nio is already installed:

if python module isn't loaded yet, then do: 
> module load python/2.7-ve0

to start using interactively in ipython:
> ipython
> import Nio


== Other Netcdf reading/writing packages ==

 * python-netcdf
 * python-scipy


To use scipy NetCDF interface, python-scipy package should be installed. 
The NetCDF reading/writing routines can be invoked like:

{{{#!python
from scipy.io import netcdf
}}}

In general it has syntax very similar to PyNio



= GDAL =

== Install on Debian or Ubuntu ==

{{{
sudo apt-get install python-gdal 
}}}

== Reading raster files ==
Sample file: [[ftp://ftp.remotesensing.org/pub/geotiff/samples/usgs/m30dem.tif | USGS DEM  ]]

{{{#!python
import gdal
from gdalconst import *


# Step 1
# open file for reading
dataset = gdal.Open( 'm30dem.tif', GA_ReadOnly )

# check what format driver is used
print 'Driver: ', dataset.GetDriver().ShortName,'/', \
          dataset.GetDriver().LongName

# Check image size
print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \
          'x',dataset.RasterCount

# check what projection it is
print 'Projection is ',dataset.GetProjection()
    
geotransform = dataset.GetGeoTransform()
if not geotransform is None:
    print 'Origin = (',geotransform[0], ',',geotransform[3],')'
    print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'

# Step 2
# Fetch raster band
band = dataset.GetRasterBand(1)
print 'Band Type=',gdal.GetDataTypeName(band.DataType)
min = band.GetMinimum()
max = band.GetMaximum()
if min is None or max is None:
    (min,max) = band.ComputeRasterMinMax(1)
    print 'Min=%.3f, Max=%.3f' % (min,max)

if band.GetOverviewCount() > 0:
    print 'Band has ', band.GetOverviewCount(), ' overviews.'

if not band.GetRasterColorTable() is None:
    print 'Band has a color table with ', \
    band.GetRasterColorTable().GetCount(), ' entries.'


# Step 3
# Read raster band contents
scanline = band.ReadRaster( 0, 0, band.XSize, 1, \
                                     band.XSize, 1, GDT_Float32 )

}}}


== Convert coordinate system == 


With Gdal and Osr it is possible to convert affine coordinates into lan/lon (example taken from [[ http://gis.stackexchange.com/questions/8392/how-do-i-convert-affine-coordinates-to-lat-lng | gis.stachexchange.com ]]


{{{#!python
import osr
import gdal

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

}}}

Same can be achieved with {{{gdal_translate}}} utility
{{{
gdal_translate -a_srs epsg:4326 srcfile new_file
 }}}



= Basemap =

{{{#!python
# Plot CCM precipitation data with basemap

import Nio
from mpl_toolkits.basemap import Basemap

f = Nio.open_file('sresa1b_ncar_ccsm3_0_run1_200001.nc', 'r')

lat,lon = f.variables['lat'][:],f.variables['lon'][:]
lat.shape

# note that lat and lon variables are 1D, we need 2D grid
lon_g,lat_g = numpy.meshgrid(lon,lat)

# next few lines from the Basemap tutorial
import matplotlib.pyplot as plt

# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
# are the lat/lon values of the lower left and upper right corners
# of the map.
# lat_ts is the latitude of true scale.
# resolution = 'c' means use crude resolution coastlines.
m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80,\
            llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')

plt.title("Mercator Projection")


# convert latitude and longitude into map coordinates
x,y = m(lon_g,lat_g)

m.pcolormesh(x,y,pr)


m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,91.,30.))
m.drawmeridians(np.arange(-180.,181.,60.))
m.drawmapboundary(fill_color='aqua')

plt.show()

}}}