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

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

Added extra file conversion classes to file_conversion.

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