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

Last change on this file since 7814 was 7814, checked in by hudson, 14 years ago

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.

File size: 4.5 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(name_in, stencil, cellsize_new, name_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    if name_in[-4:] != '.dem':
29        raise IOError('Input file %s should be of type .dem.' % name_in)
30
31    if name_out != None and basename_out[-4:] != '.dem':
32        raise IOError('Input file %s should be of type .dem.' % name_out)
33
34    #Open existing netcdf file to read
35    infile = NetCDFFile(name_in, netcdf_mode_r)
36
37    if verbose: log.critical('Reading DEM from %s' % inname)
38
39    # Read metadata (convert from numpy.int32 to int where appropriate)
40    ncols = int(infile.ncols[0])
41    nrows = int(infile.nrows[0])
42    xllcorner = infile.xllcorner[0]
43    yllcorner = infile.yllcorner[0]
44    cellsize = int(infile.cellsize[0])
45    NODATA_value = int(infile.NODATA_value[0])
46    zone = int(infile.zone[0])
47    false_easting = infile.false_easting[0]
48    false_northing = infile.false_northing[0]
49    projection = infile.projection
50    datum = infile.datum
51    units = infile.units
52
53    dem_elevation = infile.variables['elevation']
54
55    #Get output file name
56    if name_out == None:
57        outname = name_in[:-4] + '_' + repr(cellsize_new) + '.dem'
58    else:
59        outname = name_out
60
61    if verbose: log.critical('Write decimated NetCDF file to %s' % outname)
62
63    #Determine some dimensions for decimated grid
64    (nrows_stencil, ncols_stencil) = stencil.shape
65    x_offset = ncols_stencil / 2
66    y_offset = nrows_stencil / 2
67    cellsize_ratio = int(cellsize_new / cellsize)
68    ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio
69    nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio
70
71    #print type(ncols_new), ncols_new
72   
73    #Open netcdf file for output
74    outfile = NetCDFFile(outname, netcdf_mode_w)
75
76    #Create new file
77    outfile.institution = 'Geoscience Australia'
78    outfile.description = 'NetCDF DEM format for compact and portable ' \
79                          'storage of spatial point data'
80
81    #Georeferencing
82    outfile.zone = zone
83    outfile.projection = projection
84    outfile.datum = datum
85    outfile.units = units
86
87    outfile.cellsize = cellsize_new
88    outfile.NODATA_value = NODATA_value
89    outfile.false_easting = false_easting
90    outfile.false_northing = false_northing
91
92    outfile.xllcorner = xllcorner + (x_offset * cellsize)
93    outfile.yllcorner = yllcorner + (y_offset * cellsize)
94    outfile.ncols = ncols_new
95    outfile.nrows = nrows_new
96
97    # dimension definition
98    #print nrows_new, ncols_new, nrows_new*ncols_new
99    #print type(nrows_new), type(ncols_new), type(nrows_new*ncols_new)
100    outfile.createDimension('number_of_points', nrows_new*ncols_new)
101
102    # variable definition
103    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
104
105    # Get handle to the variable
106    elevation = outfile.variables['elevation']
107
108    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
109
110    #Store data
111    global_index = 0
112    for i in range(nrows_new):
113        if verbose: log.critical('Processing row %d of %d' % (i, nrows_new))
114
115        lower_index = global_index
116        telev = num.zeros(ncols_new, num.float)
117        local_index = 0
118        trow = i * cellsize_ratio
119
120        for j in range(ncols_new):
121            tcol = j * cellsize_ratio
122            tmp = dem_elevation_r[trow:trow+nrows_stencil,
123                                  tcol:tcol+ncols_stencil]
124
125            #if dem contains 1 or more NODATA_values set value in
126            #decimated dem to NODATA_value, else compute decimated
127            #value using stencil
128            if num.sum(num.sum(num.equal(tmp, NODATA_value))) > 0:
129                telev[local_index] = NODATA_value
130            else:
131                telev[local_index] = num.sum(num.sum(tmp * stencil))
132
133            global_index += 1
134            local_index += 1
135
136        upper_index = global_index
137
138        elevation[lower_index:upper_index] = telev
139
140    assert global_index == nrows_new*ncols_new, \
141           'index not equal to number of points'
142
143    infile.close()
144    outfile.close()
145
Note: See TracBrowser for help on using the repository browser.