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

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

Added extra file conversion classes to file_conversion.

File size: 5.1 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(basename_in,
41                              basename_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    root = basename_in
102
103    # Get ASCII file
104    infile = open(root + '.sdf', 'r')
105
106    if verbose: log.critical('Reading DEM from %s' % (root + '.sdf'))
107
108    lines = infile.readlines()
109    infile.close()
110
111    if verbose: log.critical('Converting to pts format')
112
113    # Scan through the header, picking up stuff we need.
114    i = 0
115    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
116        i += 1
117
118    assert lines[i].strip().upper() == 'BEGIN HEADER:'
119    i += 1
120
121    assert lines[i].strip().upper().startswith('UNITS:')
122    units = lines[i].strip().split()[1]
123    i += 1
124
125    assert lines[i].strip().upper().startswith('DTM TYPE:')
126    i += 1
127
128    assert lines[i].strip().upper().startswith('DTM:')
129    i += 1
130
131    assert lines[i].strip().upper().startswith('STREAM')
132    i += 1
133
134    assert lines[i].strip().upper().startswith('CROSS')
135    i += 1
136
137    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
138    projection = lines[i].strip().split(':')[1]
139    i += 1
140
141    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
142    zone = int(lines[i].strip().split(':')[1])
143    i += 1
144
145    assert lines[i].strip().upper().startswith('DATUM:')
146    datum = lines[i].strip().split(':')[1]
147    i += 1
148
149    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
150    i += 1
151
152    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
153    i += 1
154
155    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
156    number_of_cross_sections = int(lines[i].strip().split(':')[1])
157    i += 1
158
159    # Now read all points
160    points = []
161    elevation = []
162    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
163        for k, entry in enumerate(entries):
164            points.append(entry[:2])
165            elevation.append(entry[2])
166
167    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
168          %(j+1, number_of_cross_sections)
169    assert j+1 == number_of_cross_sections, msg
170
171    # Get output file, write PTS data
172    if basename_out == None:
173        ptsname = root + '.pts'
174    else:
175        ptsname = basename_out + '.pts'
176
177    geo_ref = Geo_reference(zone, 0, 0, datum, projection, units)
178    geo = Geospatial_data(points, {"elevation":elevation},
179                          verbose=verbose, geo_reference=geo_ref)
180    geo.export_points_file(ptsname)
Note: See TracBrowser for help on using the repository browser.