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

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

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

File size: 5.3 KB
Line 
1
2from anuga.coordinate_transforms.geo_reference import Geo_reference
3from anuga.geospatial_data.geospatial_data import Geospatial_data
4
5##
6# @brief  Return block of surface lines for each cross section
7# @param lines Iterble  of text lines to process.
8# @note BROKEN?  UNUSED?
9def _read_hecras_cross_sections(lines):
10    """Return block of surface lines for each cross section
11    Starts with SURFACE LINE,
12    Ends with END CROSS-SECTION
13    """
14
15    points = []
16
17    reading_surface = False
18    for i, line in enumerate(lines):
19        if len(line.strip()) == 0:    #Ignore blanks
20            continue
21
22        if lines[i].strip().startswith('SURFACE LINE'):
23            reading_surface = True
24            continue
25
26        if lines[i].strip().startswith('END') and reading_surface:
27            yield points
28            reading_surface = False
29            points = []
30
31        if reading_surface:
32            fields = line.strip().split(',')
33            easting = float(fields[0])
34            northing = float(fields[1])
35            elevation = float(fields[2])
36            points.append([easting, northing, elevation])
37
38
39
40def sdf2pts(name_in,
41                              name_out=None,
42                              verbose=False):
43    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
44
45# @brief Convert HECRAS (.sdf) file to PTS file.
46# @param basename_in Sterm of input filename.
47# @param basename_out Sterm of output filename.
48# @param verbose True if this function is to be verbose.
49
50    Example:
51
52# RAS export file created on Mon 15Aug2005 11:42
53# by HEC-RAS Version 3.1.1
54
55BEGIN HEADER:
56  UNITS: METRIC
57  DTM TYPE: TIN
58  DTM: v:\1\cit\perth_topo\river_tin
59  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
60  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
61  MAP PROJECTION: UTM
62  PROJECTION ZONE: 50
63  DATUM: AGD66
64  VERTICAL DATUM:
65  NUMBER OF REACHES:  19
66  NUMBER OF CROSS-SECTIONS:  14206
67END HEADER:
68
69Only the SURFACE LINE data of the following form will be utilised
70  CROSS-SECTION:
71    STREAM ID:Southern-Wungong
72    REACH ID:Southern-Wungong
73    STATION:19040.*
74    CUT LINE:
75      405548.671603161 , 6438142.7594925
76      405734.536092045 , 6438326.10404912
77      405745.130459356 , 6438331.48627354
78      405813.89633823 , 6438368.6272789
79    SURFACE LINE:
80     405548.67,   6438142.76,   35.37
81     405552.24,   6438146.28,   35.41
82     405554.78,   6438148.78,   35.44
83     405555.80,   6438149.79,   35.44
84     405559.37,   6438153.31,   35.45
85     405560.88,   6438154.81,   35.44
86     405562.93,   6438156.83,   35.42
87     405566.50,   6438160.35,   35.38
88     405566.99,   6438160.83,   35.37
89     ...
90   END CROSS-SECTION
91
92    Convert to NetCDF pts format which is
93
94    points:  (Nx2) float array
95    elevation: N float array
96    """
97
98    import os
99    from Scientific.IO.NetCDF import NetCDFFile
100
101    if name_in[-4:] != '.sdf':
102        raise IOError('Input file %s should be of type .sdf.' % name_in)
103
104    if name_out is None:
105        name_out = name_in[:-4] + '.pts'
106    elif name_out[-4:] != '.pts':
107        raise IOError('Input file %s should be of type .pts.' % name_out)
108
109    # Get ASCII file
110    infile = open(name_in, 'r')
111
112    if verbose: log.critical('Reading DEM from %s' % (root + '.sdf'))
113
114    lines = infile.readlines()
115    infile.close()
116
117    if verbose: log.critical('Converting to pts format')
118
119    # Scan through the header, picking up stuff we need.
120    i = 0
121    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
122        i += 1
123
124    assert lines[i].strip().upper() == 'BEGIN HEADER:'
125    i += 1
126
127    assert lines[i].strip().upper().startswith('UNITS:')
128    units = lines[i].strip().split()[1]
129    i += 1
130
131    assert lines[i].strip().upper().startswith('DTM TYPE:')
132    i += 1
133
134    assert lines[i].strip().upper().startswith('DTM:')
135    i += 1
136
137    assert lines[i].strip().upper().startswith('STREAM')
138    i += 1
139
140    assert lines[i].strip().upper().startswith('CROSS')
141    i += 1
142
143    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
144    projection = lines[i].strip().split(':')[1]
145    i += 1
146
147    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
148    zone = int(lines[i].strip().split(':')[1])
149    i += 1
150
151    assert lines[i].strip().upper().startswith('DATUM:')
152    datum = lines[i].strip().split(':')[1]
153    i += 1
154
155    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
156    i += 1
157
158    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
159    i += 1
160
161    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
162    number_of_cross_sections = int(lines[i].strip().split(':')[1])
163    i += 1
164
165    # Now read all points
166    points = []
167    elevation = []
168    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
169        for k, entry in enumerate(entries):
170            points.append(entry[:2])
171            elevation.append(entry[2])
172
173    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
174          %(j+1, number_of_cross_sections)
175    assert j+1 == number_of_cross_sections, msg
176
177    # Get output file, write PTS data
178    if name_out == None:
179        ptsname = name_in[:-4] + '.pts'
180    else:
181        ptsname = name_out
182
183    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
184    geo = Geospatial_data(points, {"elevation":elevation},
185                          verbose=verbose, geo_reference=geo_ref)
186    geo.export_points_file(ptsname)
Note: See TracBrowser for help on using the repository browser.