1 # load Nio module
   2 import Nio
   3 
   4 # open netcdf file
   5 f = Nio.open_file('sresa1b_ncar_ccsm3_0_run1_200001.nc', 'r') # 'r' stands for "read rights"
   6 
   7 # check contents
   8 print f
   9 
  10 # get dimensions
  11 dimNames = f.dimensions.keys()
  12 
  13 # get the size of dimension
  14 dimSize = f.dimensions['dimName']
  15 
  16 # read single variable:
  17 pr = f.variables['pr']
  18 print pr
  19 
  20 # read variables contents
  21 pr_data = pr[:]
  22 print pr_data
  23 
  24 # print slice of the variable data
  25 slice = pr_data[0,::-1,:]
  26 
  27 # close file
  28 f.close()

   1 # create new file
   2 import Nio
   3 import numpy
   4 
   5 f = Nio.open_file('test.nc', 'w') # "w" stands for writing rights
   6 
   7 # create a dimension
   8 f.create_dimension('time',100)
   9 
  10 # create variable
  11 f.create_variable('temperature', 'i', ('time',)) # dimension 'time' must be created prior to this step
  12 
  13 # dimensions and variables can have same names
  14 f.create_variable('time', 'i', ('time',))
  15 
  16 # Dimension types:
  17 #    'd': 64 bit float
  18 #    'f': 32 bit float
  19 #    'l': long
  20 #    'i': 32 bit integer
  21 #    'h': 16 bit integer
  22 #    'b': 8 bit integer
  23 #    'S1': character 
  24 
  25 # Create an attribute 
  26 f.variables['temperature'].units = "K"
  27 
  28 # Add scaling factor and an offset
  29 # when the variable will be read, it will need to be multiplied by the scale_factor first and the added to the offset value
  30 f.variables['temperature'].scale_factor = 0.1
  31 f.variables['temperature'].add_offset = 273
  32 
  33 # create variable contents and assign it to the NetCDF variable
  34 temp = numpy.arange(0,100,1, dtype = numpy.int32)
  35 f.variables['temperature'][:] = temp
  36 
  37 time = numpy.arange(1000,1100,1, dtype = numpy.int32)
  38 f.variables['time'][:] = time
  39 
  40 f.close()

   1 from scipy.io import netcdf

   1 import gdal
   2 from gdalconst import *
   3 
   4 
   5 # Step 1
   6 # open file for reading
   7 dataset = gdal.Open( 'm30dem.tif', GA_ReadOnly )
   8 
   9 # check what format driver is used
  10 print 'Driver: ', dataset.GetDriver().ShortName,'/', \
  11           dataset.GetDriver().LongName
  12 
  13 # Check image size
  14 print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \
  15           'x',dataset.RasterCount
  16 
  17 # check what projection it is
  18 print 'Projection is ',dataset.GetProjection()
  19     
  20 geotransform = dataset.GetGeoTransform()
  21 if not geotransform is None:
  22     print 'Origin = (',geotransform[0], ',',geotransform[3],')'
  23     print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'
  24 
  25 # Step 2
  26 # Fetch raster band
  27 band = dataset.GetRasterBand(1)
  28 print 'Band Type=',gdal.GetDataTypeName(band.DataType)
  29 min = band.GetMinimum()
  30 max = band.GetMaximum()
  31 if min is None or max is None:
  32     (min,max) = band.ComputeRasterMinMax(1)
  33     print 'Min=%.3f, Max=%.3f' % (min,max)
  34 
  35 if band.GetOverviewCount() > 0:
  36     print 'Band has ', band.GetOverviewCount(), ' overviews.'
  37 
  38 if not band.GetRasterColorTable() is None:
  39     print 'Band has a color table with ', \
  40     band.GetRasterColorTable().GetCount(), ' entries.'
  41 
  42 
  43 # Step 3
  44 # Read raster band contents
  45 scanline = band.ReadRaster( 0, 0, band.XSize, 1, \
  46                                      band.XSize, 1, GDT_Float32 )

   1 import osr
   2 import gdal
   3 
   4 srs = osr.SpatialReference()
   5 srs.ImportFromWkt(dataset.GetProjection())
   6 
   7 srsLatLong = srs.CloneGeogCS()
   8 ct = osr.CoordinateTransformation(srs,srsLatLong)
   9 print ct.TransformPoint(originX,originY)

   1 # Plot CCM precipitation data with basemap
   2 
   3 import Nio
   4 from mpl_toolkits.basemap import Basemap
   5 
   6 f = Nio.open_file('sresa1b_ncar_ccsm3_0_run1_200001.nc', 'r')
   7 
   8 lat,lon = f.variables['lat'][:],f.variables['lon'][:]
   9 lat.shape
  10 
  11 # note that lat and lon variables are 1D, we need 2D grid
  12 lon_g,lat_g = numpy.meshgrid(lon,lat)
  13 
  14 # next few lines from the Basemap tutorial
  15 import matplotlib.pyplot as plt
  16 
  17 # llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
  18 # are the lat/lon values of the lower left and upper right corners
  19 # of the map.
  20 # lat_ts is the latitude of true scale.
  21 # resolution = 'c' means use crude resolution coastlines.
  22 m = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80,\
  23             llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
  24 
  25 plt.title("Mercator Projection")
  26 
  27 
  28 # convert latitude and longitude into map coordinates
  29 x,y = m(lon_g,lat_g)
  30 
  31 m.pcolormesh(x,y,pr)
  32 
  33 
  34 m.drawcoastlines()
  35 m.fillcontinents(color='coral',lake_color='aqua')
  36 # draw parallels and meridians.
  37 m.drawparallels(np.arange(-90.,91.,30.))
  38 m.drawmeridians(np.arange(-180.,181.,60.))
  39 m.drawmapboundary(fill_color='aqua')
  40 
  41 plt.show()