source: trunk/anuga_core/source/anuga/file_conversion/asc2dem.py @ 7804

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

Fixed up failing tests, updated user guide with new API (first few chapters only).
File size: 6.7 KB
Line 
1# external modules
2import numpy as num
3import anuga.utilities.log as log
4
5# ANUGA modules
6from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
7                            netcdf_float
8
9##
10# @brief Convert ASC file to DEM file.
11# formerly convert_dem_from_ascii2netcdf
12# @param basename_in Stem of input filename.
13# @param basename_out Stem of output filename.
14# @param use_cache ??
15# @param verbose True if this function is to be verbose.
16# @return
17# @note A PRJ file with same stem basename must exist and is used to fix the
18#       UTM zone, datum, false northings and eastings.
19def asc2dem(basename_in, basename_out=None,
20                                  use_cache=False,
21                                  verbose=False):
22    """Read Digital Elevation model from the following ASCII format (.asc)
23
24    Example:
25    ncols         3121
26    nrows         1800
27    xllcorner     722000
28    yllcorner     5893000
29    cellsize      25
30    NODATA_value  -9999
31    138.3698 137.4194 136.5062 135.5558 ..........
32
33    Convert basename_in + '.asc' to NetCDF format (.dem)
34    mimicking the ASCII format closely.
35
36    An accompanying file with same basename_in but extension .prj must exist
37    and is used to fix the UTM zone, datum, false northings and eastings.
38
39    The prj format is assumed to be as
40
41    Projection    UTM
42    Zone          56
43    Datum         WGS84
44    Zunits        NO
45    Units         METERS
46    Spheroid      WGS84
47    Xshift        0.0000000000
48    Yshift        10000000.0000000000
49    Parameters
50    """
51
52    kwargs = {'basename_out': basename_out, 'verbose': verbose}
53
54    if use_cache is True:
55        from caching import cache
56        result = cache(_convert_dem_from_ascii2netcdf, basename_in, kwargs,
57                       dependencies=[basename_in + '.asc',
58                                     basename_in + '.prj'],
59                       verbose=verbose)
60
61    else:
62        result = apply(_convert_dem_from_ascii2netcdf, [basename_in], kwargs)
63
64    return result
65
66
67##
68# @brief Convert an ASC file to a DEM file.
69# @param basename_in Stem of input filename.
70# @param basename_out Stem of output filename.
71# @param verbose True if this function is to be verbose.
72def _convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
73                                   verbose = False):
74    """Read Digital Elevation model from the following ASCII format (.asc)
75
76    Internal function. See public function convert_dem_from_ascii2netcdf
77    for details.
78    """
79
80    import os
81    from Scientific.IO.NetCDF import NetCDFFile
82
83    root = basename_in
84
85    # Read Meta data
86    if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
87
88    metadatafile = open(root + '.prj')
89    metalines = metadatafile.readlines()
90    metadatafile.close()
91
92    L = metalines[0].strip().split()
93    assert L[0].strip().lower() == 'projection'
94    projection = L[1].strip()                   #TEXT
95
96    L = metalines[1].strip().split()
97    assert L[0].strip().lower() == 'zone'
98    zone = int(L[1].strip())
99
100    L = metalines[2].strip().split()
101    assert L[0].strip().lower() == 'datum'
102    datum = L[1].strip()                        #TEXT
103
104    L = metalines[3].strip().split()
105    assert L[0].strip().lower() == 'zunits'     #IGNORE
106    zunits = L[1].strip()                       #TEXT
107
108    L = metalines[4].strip().split()
109    assert L[0].strip().lower() == 'units'
110    units = L[1].strip()                        #TEXT
111
112    L = metalines[5].strip().split()
113    assert L[0].strip().lower() == 'spheroid'   #IGNORE
114    spheroid = L[1].strip()                     #TEXT
115
116    L = metalines[6].strip().split()
117    assert L[0].strip().lower() == 'xshift'
118    false_easting = float(L[1].strip())
119
120    L = metalines[7].strip().split()
121    assert L[0].strip().lower() == 'yshift'
122    false_northing = float(L[1].strip())
123
124    #Read DEM data
125    datafile = open(basename_in + '.asc')
126
127    if verbose: log.critical('Reading DEM from %s' % (basename_in + '.asc'))
128
129    lines = datafile.readlines()
130    datafile.close()
131
132    if verbose: log.critical('Got %d lines' % len(lines))
133
134    ncols = int(lines[0].split()[1].strip())
135    nrows = int(lines[1].split()[1].strip())
136
137    # Do cellsize (line 4) before line 2 and 3
138    cellsize = float(lines[4].split()[1].strip())
139
140    # Checks suggested by Joaquim Luis
141    # Our internal representation of xllcorner
142    # and yllcorner is non-standard.
143    xref = lines[2].split()
144    if xref[0].strip() == 'xllcorner':
145        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
146    elif xref[0].strip() == 'xllcenter':
147        xllcorner = float(xref[1].strip())
148    else:
149        msg = 'Unknown keyword: %s' % xref[0].strip()
150        raise Exception, msg
151
152    yref = lines[3].split()
153    if yref[0].strip() == 'yllcorner':
154        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
155    elif yref[0].strip() == 'yllcenter':
156        yllcorner = float(yref[1].strip())
157    else:
158        msg = 'Unknown keyword: %s' % yref[0].strip()
159        raise Exception, msg
160
161    NODATA_value = int(lines[5].split()[1].strip())
162
163    assert len(lines) == nrows + 6
164
165    if basename_out == None:
166        netcdfname = root + '.dem'
167    else:
168        netcdfname = basename_out + '.dem'
169
170    if verbose: log.critical('Store to NetCDF file %s' % netcdfname)
171
172    # NetCDF file definition
173    fid = NetCDFFile(netcdfname, netcdf_mode_w)
174
175    #Create new file
176    fid.institution = 'Geoscience Australia'
177    fid.description = 'NetCDF DEM format for compact and portable storage ' \
178                      'of spatial point data'
179
180    fid.ncols = ncols
181    fid.nrows = nrows
182    fid.xllcorner = xllcorner
183    fid.yllcorner = yllcorner
184    fid.cellsize = cellsize
185    fid.NODATA_value = NODATA_value
186
187    fid.zone = zone
188    fid.false_easting = false_easting
189    fid.false_northing = false_northing
190    fid.projection = projection
191    fid.datum = datum
192    fid.units = units
193
194    # dimension definitions
195    fid.createDimension('number_of_rows', nrows)
196    fid.createDimension('number_of_columns', ncols)
197
198    # variable definitions
199    fid.createVariable('elevation', netcdf_float, ('number_of_rows',
200                                                   'number_of_columns'))
201
202    # Get handles to the variables
203    elevation = fid.variables['elevation']
204
205    #Store data
206    n = len(lines[6:])
207    for i, line in enumerate(lines[6:]):
208        fields = line.split()
209        if verbose and i % ((n+10)/10) == 0:
210            log.critical('Processing row %d of %d' % (i, nrows))
211           
212        if len(fields) != ncols:
213            msg = 'Wrong number of columns in file "%s" line %d\n' % (basename_in + '.asc', i)
214            msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
215            raise Exception, msg
216
217        elevation[i, :] = num.array([float(x) for x in fields])
218
219    fid.close()
Note: See TracBrowser for help on using the repository browser.