Size: 4518
Comment:
|
Size: 5236
Comment:
|
Deletions are marked like this. | Additions are marked like this. |
Line 140: | Line 140: |
# Step 1 |
|
Line 159: | Line 161: |
}}} | # 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 ) }}} |
Line 162: | Line 188: |
Contents
PyNIO
Find full PyNio documentation here
Download a 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
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()
Create NetCDF file
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()
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
GDAL
Install on Debian or Ubuntu
sudo apt-get install python-gdal
Reading raster files
e.g. a ftp://ftp.remotesensing.org/pub/geotiff/samples/usgs/m30dem.tif
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 )
== 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
Same can be achieved with gdal_translate utility
gdal_translate -a_srs epsg:4326 srcfile new_file