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

Last change on this file since 8150 was 8150, checked in by wilsonr, 13 years ago

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