source: trunk/anuga_core/source/anuga/file_conversion/dem2dem.py @ 7778

Last change on this file since 7778 was 7778, checked in by hudson, 13 years ago

Cleaned up unit tests to use API.

File size: 4.3 KB
Line 
1""" Load a DEM file, decimate it, and resave it.
2"""
3
4import numpy as num
5
6from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_float
7
8def dem2dem(basename_in, stencil, cellsize_new, basename_out=None,
9                 verbose=False):
10    """Read Digitial Elevation model from the following NetCDF format (.dem)
11
12    Example:
13
14    ncols         3121
15    nrows         1800
16    xllcorner     722000
17    yllcorner     5893000
18    cellsize      25
19    NODATA_value  -9999
20    138.3698 137.4194 136.5062 135.5558 ..........
21
22    Decimate data to cellsize_new using stencil and write to NetCDF dem format.
23    """
24
25    import os
26    from Scientific.IO.NetCDF import NetCDFFile
27
28    root = basename_in
29    inname = root + '.dem'
30
31    #Open existing netcdf file to read
32    infile = NetCDFFile(inname, netcdf_mode_r)
33
34    if verbose: log.critical('Reading DEM from %s' % inname)
35
36    # Read metadata (convert from numpy.int32 to int where appropriate)
37    ncols = int(infile.ncols[0])
38    nrows = int(infile.nrows[0])
39    xllcorner = infile.xllcorner[0]
40    yllcorner = infile.yllcorner[0]
41    cellsize = int(infile.cellsize[0])
42    NODATA_value = int(infile.NODATA_value[0])
43    zone = int(infile.zone[0])
44    false_easting = infile.false_easting[0]
45    false_northing = infile.false_northing[0]
46    projection = infile.projection
47    datum = infile.datum
48    units = infile.units
49
50    dem_elevation = infile.variables['elevation']
51
52    #Get output file name
53    if basename_out == None:
54        outname = root + '_' + repr(cellsize_new) + '.dem'
55    else:
56        outname = basename_out + '.dem'
57
58    if verbose: log.critical('Write decimated NetCDF file to %s' % outname)
59
60    #Determine some dimensions for decimated grid
61    (nrows_stencil, ncols_stencil) = stencil.shape
62    x_offset = ncols_stencil / 2
63    y_offset = nrows_stencil / 2
64    cellsize_ratio = int(cellsize_new / cellsize)
65    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
66    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
67
68    #print type(ncols_new), ncols_new
69   
70    #Open netcdf file for output
71    outfile = NetCDFFile(outname, netcdf_mode_w)
72
73    #Create new file
74    outfile.institution = 'Geoscience Australia'
75    outfile.description = 'NetCDF DEM format for compact and portable ' \
76                          'storage of spatial point data'
77
78    #Georeferencing
79    outfile.zone = zone
80    outfile.projection = projection
81    outfile.datum = datum
82    outfile.units = units
83
84    outfile.cellsize = cellsize_new
85    outfile.NODATA_value = NODATA_value
86    outfile.false_easting = false_easting
87    outfile.false_northing = false_northing
88
89    outfile.xllcorner = xllcorner + (x_offset * cellsize)
90    outfile.yllcorner = yllcorner + (y_offset * cellsize)
91    outfile.ncols = ncols_new
92    outfile.nrows = nrows_new
93
94    # dimension definition
95    #print nrows_new, ncols_new, nrows_new*ncols_new
96    #print type(nrows_new), type(ncols_new), type(nrows_new*ncols_new)
97    outfile.createDimension('number_of_points', nrows_new*ncols_new)
98
99    # variable definition
100    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
101
102    # Get handle to the variable
103    elevation = outfile.variables['elevation']
104
105    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
106
107    #Store data
108    global_index = 0
109    for i in range(nrows_new):
110        if verbose: log.critical('Processing row %d of %d' % (i, nrows_new))
111
112        lower_index = global_index
113        telev = num.zeros(ncols_new, num.float)
114        local_index = 0
115        trow = i * cellsize_ratio
116
117        for j in range(ncols_new):
118            tcol = j * cellsize_ratio
119            tmp = dem_elevation_r[trow:trow+nrows_stencil,
120                                  tcol:tcol+ncols_stencil]
121
122            #if dem contains 1 or more NODATA_values set value in
123            #decimated dem to NODATA_value, else compute decimated
124            #value using stencil
125            if num.sum(num.sum(num.equal(tmp, NODATA_value))) > 0:
126                telev[local_index] = NODATA_value
127            else:
128                telev[local_index] = num.sum(num.sum(tmp * stencil))
129
130            global_index += 1
131            local_index += 1
132
133        upper_index = global_index
134
135        elevation[lower_index:upper_index] = telev
136
137    assert global_index == nrows_new*ncols_new, \
138           'index not equal to number of points'
139
140    infile.close()
141    outfile.close()
142
Note: See TracBrowser for help on using the repository browser.