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

Last change on this file since 7866 was 7814, checked in by James Hudson, 15 years ago

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.
File size: 6.4 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
57##
58# @brief Convert an ASC file to a DEM file.
59# @param basename_in Stem of input filename.
60# @param basename_out Stem of output filename.
61# @param verbose True if this function is to be verbose.
62def _convert_dem_from_ascii2netcdf(name_in, name_out = None,
63                                   verbose = False):
64    """Read Digital Elevation model from the following ASCII format (.asc)
65
66    Internal function. See public function convert_dem_from_ascii2netcdf
67    for details.
68    """
69
70    import os
71    from Scientific.IO.NetCDF import NetCDFFile
72
73    root = name_in[:-4]
74
75    # Read Meta data
76    if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
77
78    metadatafile = open(root + '.prj')
79    metalines = metadatafile.readlines()
80    metadatafile.close()
81
82    L = metalines[0].strip().split()
83    assert L[0].strip().lower() == 'projection'
84    projection = L[1].strip()                   #TEXT
85
86    L = metalines[1].strip().split()
87    assert L[0].strip().lower() == 'zone'
88    zone = int(L[1].strip())
89
90    L = metalines[2].strip().split()
91    assert L[0].strip().lower() == 'datum'
92    datum = L[1].strip()                        #TEXT
93
94    L = metalines[3].strip().split()
95    assert L[0].strip().lower() == 'zunits'     #IGNORE
96    zunits = L[1].strip()                       #TEXT
97
98    L = metalines[4].strip().split()
99    assert L[0].strip().lower() == 'units'
100    units = L[1].strip()                        #TEXT
101
102    L = metalines[5].strip().split()
103    assert L[0].strip().lower() == 'spheroid'   #IGNORE
104    spheroid = L[1].strip()                     #TEXT
105
106    L = metalines[6].strip().split()
107    assert L[0].strip().lower() == 'xshift'
108    false_easting = float(L[1].strip())
109
110    L = metalines[7].strip().split()
111    assert L[0].strip().lower() == 'yshift'
112    false_northing = float(L[1].strip())
113
114    if name_in[-4:] != '.asc':
115        raise IOError('Input file %s should be of type .asc.' % name_in)
116
117    #Read DEM data
118    datafile = open(name_in)
119
120    if verbose: log.critical('Reading DEM from %s' % (name_in))
121
122    lines = datafile.readlines()
123    datafile.close()
124
125    if verbose: log.critical('Got %d lines' % len(lines))
126
127    ncols = int(lines[0].split()[1].strip())
128    nrows = int(lines[1].split()[1].strip())
129
130    # Do cellsize (line 4) before line 2 and 3
131    cellsize = float(lines[4].split()[1].strip())
132
133    # Checks suggested by Joaquim Luis
134    # Our internal representation of xllcorner
135    # and yllcorner is non-standard.
136    xref = lines[2].split()
137    if xref[0].strip() == 'xllcorner':
138        xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
139    elif xref[0].strip() == 'xllcenter':
140        xllcorner = float(xref[1].strip())
141    else:
142        msg = 'Unknown keyword: %s' % xref[0].strip()
143        raise Exception, msg
144
145    yref = lines[3].split()
146    if yref[0].strip() == 'yllcorner':
147        yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
148    elif yref[0].strip() == 'yllcenter':
149        yllcorner = float(yref[1].strip())
150    else:
151        msg = 'Unknown keyword: %s' % yref[0].strip()
152        raise Exception, msg
153
154    NODATA_value = int(lines[5].split()[1].strip())
155
156    assert len(lines) == nrows + 6
157
158    if name_out == None:
159        netcdfname = name_in[:-4]+'.dem'
160    else:
161        netcdfname = name_out + '.dem'
162
163    if verbose: log.critical('Store to NetCDF file %s' % netcdfname)
164
165    # NetCDF file definition
166    fid = NetCDFFile(netcdfname, netcdf_mode_w)
167
168    #Create new file
169    fid.institution = 'Geoscience Australia'
170    fid.description = 'NetCDF DEM format for compact and portable storage ' \
171                      'of spatial point data'
172
173    fid.ncols = ncols
174    fid.nrows = nrows
175    fid.xllcorner = xllcorner
176    fid.yllcorner = yllcorner
177    fid.cellsize = cellsize
178    fid.NODATA_value = NODATA_value
179
180    fid.zone = zone
181    fid.false_easting = false_easting
182    fid.false_northing = false_northing
183    fid.projection = projection
184    fid.datum = datum
185    fid.units = units
186
187    # dimension definitions
188    fid.createDimension('number_of_rows', nrows)
189    fid.createDimension('number_of_columns', ncols)
190
191    # variable definitions
192    fid.createVariable('elevation', netcdf_float, ('number_of_rows',
193                                                   'number_of_columns'))
194
195    # Get handles to the variables
196    elevation = fid.variables['elevation']
197
198    #Store data
199    n = len(lines[6:])
200    for i, line in enumerate(lines[6:]):
201        fields = line.split()
202        if verbose and i % ((n+10)/10) == 0:
203            log.critical('Processing row %d of %d' % (i, nrows))
204           
205        if len(fields) != ncols:
206            msg = 'Wrong number of columns in file "%s" line %d\n' % (name_in, i)
207            msg += 'I got %d elements, but there should have been %d\n' % (len(fields), ncols)
208            raise Exception, msg
209
210        elevation[i, :] = num.array([float(x) for x in fields])
211
212    fid.close()
Note: See TracBrowser for help on using the repository browser.