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

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

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

File size: 7.8 KB
Line 
1# external modules
2import numpy as num
3
4# ANUGA modules
5import anuga.utilities.log as log
6from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \
7                            netcdf_float
8                           
9##
10# @brief Convert DEM data  to PTS data.
11# @param basename_in Stem of input filename.
12# @param basename_out Stem of output filename.
13# @param easting_min
14# @param easting_max
15# @param northing_min
16# @param northing_max
17# @param use_cache
18# @param verbose
19# @return
20def dem2pts(name_in, name_out=None,
21            easting_min=None, easting_max=None,
22            northing_min=None, northing_max=None,
23            use_cache=False, verbose=False,):
24    """Read Digitial Elevation model from the following NetCDF format (.dem)
25
26    Example:
27
28    ncols         3121
29    nrows         1800
30    xllcorner     722000
31    yllcorner     5893000
32    cellsize      25
33    NODATA_value  -9999
34    138.3698 137.4194 136.5062 135.5558 ..........
35
36    Convert to NetCDF pts format which is
37
38    points:  (Nx2) float array
39    elevation: N float array
40    """
41
42    kwargs = {'name_out': name_out,
43              'easting_min': easting_min,
44              'easting_max': easting_max,
45              'northing_min': northing_min,
46              'northing_max': northing_max,
47              'verbose': verbose}
48
49    if use_cache is True:
50        from caching import cache
51        result = cache(_dem2pts, name_in, kwargs,
52                       dependencies = [name_in],
53                       verbose = verbose)
54
55    else:
56        result = apply(_dem2pts, [name_in], kwargs)
57
58    return result
59
60
61##
62# @brief
63# @param basename_in
64# @param basename_out
65# @param verbose
66# @param easting_min
67# @param easting_max
68# @param northing_min
69# @param northing_max
70def _dem2pts(name_in, name_out=None, verbose=False,
71            easting_min=None, easting_max=None,
72            northing_min=None, northing_max=None):
73    """Read Digitial Elevation model from the following NetCDF format (.dem)
74
75    Internal function. See public function dem2pts for details.
76    """
77
78    # FIXME: Can this be written feasibly using write_pts?
79
80    import os
81    from Scientific.IO.NetCDF import NetCDFFile
82
83    if name_in[-4:] != '.dem':
84        raise IOError('Input file %s should be of type .dem.' % name_in)
85
86    if name_out != None and basename_out[-4:] != '.pts':
87        raise IOError('Input file %s should be of type .pts.' % name_out)
88
89    root = name_in[:-4]
90
91    # Get NetCDF
92    infile = NetCDFFile(name_in, netcdf_mode_r) 
93
94    if verbose: log.critical('Reading DEM from %s' % (name_in))
95
96    ncols = infile.ncols[0]
97    nrows = infile.nrows[0]
98    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
99    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
100    cellsize = infile.cellsize[0]
101    NODATA_value = infile.NODATA_value[0]
102    dem_elevation = infile.variables['elevation']
103
104    zone = infile.zone[0]
105    false_easting = infile.false_easting[0]
106    false_northing = infile.false_northing[0]
107
108    # Text strings
109    projection = infile.projection
110    datum = infile.datum
111    units = infile.units
112
113    # Get output file
114    if name_out == None:
115        ptsname = root + '.pts'
116    else:
117        ptsname = name_out
118
119    if verbose: log.critical('Store to NetCDF file %s' % ptsname)
120
121    # NetCDF file definition
122    outfile = NetCDFFile(ptsname, netcdf_mode_w)
123
124    # Create new file
125    outfile.institution = 'Geoscience Australia'
126    outfile.description = 'NetCDF pts format for compact and portable ' \
127                          'storage of spatial point data'
128
129    # Assign default values
130    if easting_min is None: easting_min = xllcorner
131    if easting_max is None: easting_max = xllcorner + ncols*cellsize
132    if northing_min is None: northing_min = yllcorner
133    if northing_max is None: northing_max = yllcorner + nrows*cellsize
134
135    # Compute offsets to update georeferencing
136    easting_offset = xllcorner - easting_min
137    northing_offset = yllcorner - northing_min
138
139    # Georeferencing
140    outfile.zone = zone
141    outfile.xllcorner = easting_min # Easting of lower left corner
142    outfile.yllcorner = northing_min # Northing of lower left corner
143    outfile.false_easting = false_easting
144    outfile.false_northing = false_northing
145
146    outfile.projection = projection
147    outfile.datum = datum
148    outfile.units = units
149
150    # Grid info (FIXME: probably not going to be used, but heck)
151    outfile.ncols = ncols
152    outfile.nrows = nrows
153
154    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
155    totalnopoints = nrows*ncols
156
157    # Calculating number of NODATA_values for each row in clipped region
158    # FIXME: use array operations to do faster
159    nn = 0
160    k = 0
161    i1_0 = 0
162    j1_0 = 0
163    thisj = 0
164    thisi = 0
165    for i in range(nrows):
166        y = (nrows-i-1)*cellsize + yllcorner
167        for j in range(ncols):
168            x = j*cellsize + xllcorner
169            if easting_min <= x <= easting_max \
170               and northing_min <= y <= northing_max:
171                thisj = j
172                thisi = i
173                if dem_elevation_r[i,j] == NODATA_value:
174                    nn += 1
175
176                if k == 0:
177                    i1_0 = i
178                    j1_0 = j
179
180                k += 1
181
182    index1 = j1_0
183    index2 = thisj
184
185    # Dimension definitions
186    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
187    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
188
189    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
190    nopoints = clippednopoints-nn
191
192    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
193
194    if verbose:
195        log.critical('There are %d values in the elevation' % totalnopoints)
196        log.critical('There are %d values in the clipped elevation'
197                     % clippednopoints)
198        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
199
200    outfile.createDimension('number_of_points', nopoints)
201    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
202
203    # Variable definitions
204    outfile.createVariable('points', netcdf_float, ('number_of_points',
205                                                    'number_of_dimensions'))
206    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
207
208    # Get handles to the variables
209    points = outfile.variables['points']
210    elevation = outfile.variables['elevation']
211
212    lenv = index2-index1+1
213
214    # Store data
215    global_index = 0
216    # for i in range(nrows):
217    for i in range(i1_0, thisi+1, 1):
218        if verbose and i % ((nrows+10)/10) == 0:
219            log.critical('Processing row %d of %d' % (i, nrows))
220
221        lower_index = global_index
222
223        v = dem_elevation_r[i,index1:index2+1]
224        no_NODATA = num.sum(v == NODATA_value)
225        if no_NODATA > 0:
226            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
227        else:
228            newcols = lenv              # ncols_in_bounding_box
229
230        telev = num.zeros(newcols, num.float)
231        tpoints = num.zeros((newcols, 2), num.float)
232
233        local_index = 0
234
235        y = (nrows-i-1)*cellsize + yllcorner
236        #for j in range(ncols):
237        for j in range(j1_0,index2+1,1):
238            x = j*cellsize + xllcorner
239            if easting_min <= x <= easting_max \
240               and northing_min <= y <= northing_max \
241               and dem_elevation_r[i,j] != NODATA_value:
242                tpoints[local_index, :] = [x-easting_min, y-northing_min]
243                telev[local_index] = dem_elevation_r[i, j]
244                global_index += 1
245                local_index += 1
246
247        upper_index = global_index
248
249        if upper_index == lower_index + newcols:
250            points[lower_index:upper_index, :] = tpoints
251            elevation[lower_index:upper_index] = telev
252
253    assert global_index == nopoints, 'index not equal to number of points'
254
255    infile.close()
256    outfile.close()
257
Note: See TracBrowser for help on using the repository browser.