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

Last change on this file since 9582 was 8832, checked in by steve, 12 years ago

Fixed some problems with dem2pts which turned up in merewether case study
File size: 10.6 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
9from asc2dem import asc2dem
10                           
11
12def dem2pts(name_in, name_out=None,
13            easting_min=None, easting_max=None,
14            northing_min=None, northing_max=None,
15            use_cache=False, verbose=False,):
16    """Read Digitial Elevation model from the following NetCDF format (.dem)
17
18    Example:
19
20    ncols         3121
21    nrows         1800
22    xllcorner     722000
23    yllcorner     5893000
24    cellsize      25
25    NODATA_value  -9999
26    138.3698 137.4194 136.5062 135.5558 ..........
27
28    name_in may be a .asc or .dem file to be converted.
29
30    Convert to NetCDF pts format which is
31
32    points:  (Nx2) float array
33    elevation: N float array
34    """
35
36    kwargs = {'name_out': name_out,
37              'easting_min': easting_min,
38              'easting_max': easting_max,
39              'northing_min': northing_min,
40              'northing_max': northing_max,
41              'verbose': verbose}
42
43    if use_cache is True:
44        from caching import cache
45        result = cache(_dem2pts, name_in, kwargs,
46                       dependencies = [name_in],
47                       verbose = verbose)
48
49    else:
50        result = apply(_dem2pts, [name_in], kwargs)
51
52    return result
53
54
55def _dem2pts(name_in, name_out=None, verbose=False,
56            easting_min=None, easting_max=None,
57            northing_min=None, northing_max=None):
58    """Read Digitial Elevation model from the following NetCDF format (.dem)
59
60    Internal function. See public function dem2pts for details.
61    """
62
63    # FIXME: Can this be written feasibly using write_pts?
64
65    import os
66    from anuga.file.netcdf import NetCDFFile
67
68    root = name_in[:-4]
69
70    if name_in[-4:] == '.asc':
71        intermediate = root + '.dem'
72        if verbose:
73            log.critical('Preconvert %s from asc to %s' % \
74                                    (name_in, intermediate))
75        asc2dem(name_in)
76        name_in = intermediate
77    elif name_in[-4:] != '.dem':
78        raise IOError('Input file %s should be of type .asc or .dem.' % name_in)
79
80    if name_out != None and basename_out[-4:] != '.pts':
81        raise IOError('Input file %s should be of type .pts.' % name_out)
82
83    # Get NetCDF
84    infile = NetCDFFile(name_in, netcdf_mode_r) 
85
86    if verbose: log.critical('Reading DEM from %s' % (name_in))
87
88    ncols = int(infile.ncols)
89    nrows = int(infile.nrows)
90    xllcorner = float(infile.xllcorner)  # Easting of lower left corner
91    yllcorner = float(infile.yllcorner)  # Northing of lower left corner
92    cellsize = float(infile.cellsize)
93    NODATA_value = float(infile.NODATA_value)
94
95    dem_elevation = infile.variables['elevation']
96
97    zone = int(infile.zone)
98    false_easting = float(infile.false_easting)
99    false_northing = float(infile.false_northing)
100
101    #print ncols, nrows, xllcorner,yllcorner, cellsize, NODATA_value, zone
102
103
104    # Text strings
105    projection = infile.projection
106    datum = infile.datum
107    units = infile.units
108
109    #print projection, datum, units
110
111    # Get output file
112    if name_out == None:
113        ptsname = root + '.pts'
114    else:
115        ptsname = name_out
116
117    if verbose: log.critical('Store to NetCDF file %s' % ptsname)
118
119    # NetCDF file definition
120    outfile = NetCDFFile(ptsname, netcdf_mode_w)
121
122    # Create new file
123    outfile.institution = 'Geoscience Australia'
124    outfile.description = 'NetCDF pts format for compact and portable ' \
125                          'storage of spatial point data'
126
127    # Assign default values
128    if easting_min is None: easting_min = xllcorner
129    if easting_max is None: easting_max = xllcorner + ncols*cellsize
130    if northing_min is None: northing_min = yllcorner
131    if northing_max is None: northing_max = yllcorner + nrows*cellsize
132
133
134    #print easting_min, easting_max, northing_min, northing_max
135
136    # Compute offsets to update georeferencing
137    easting_offset = xllcorner - easting_min
138    northing_offset = yllcorner - northing_min
139
140    # Georeferencing
141    outfile.zone = zone
142    outfile.xllcorner = easting_min # Easting of lower left corner
143    outfile.yllcorner = northing_min # Northing of lower left corner
144    outfile.false_easting = false_easting
145    outfile.false_northing = false_northing
146
147    outfile.projection = projection
148    outfile.datum = datum
149    outfile.units = units
150
151    # Grid info (FIXME: probably not going to be used, but heck)
152    outfile.ncols = ncols
153    outfile.nrows = nrows
154
155    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
156    totalnopoints = nrows*ncols
157
158
159
160
161
162
163#    #=======================================================================
164#    # Calculating number of NODATA_values for each row in clipped region
165#    # FIXME: use array operations to do faster
166#    nn = 0
167#    k = 0
168#    i1_0 = 0
169#    j1_0 = 0
170#    thisj = 0
171#    thisi = 0
172#    for i in range(nrows):
173#        y = (nrows-i-1)*cellsize + yllcorner
174#        for j in range(ncols):
175#            x = j*cellsize + xllcorner
176#            if easting_min <= x <= easting_max \
177#               and northing_min <= y <= northing_max:
178#                thisj = j
179#                thisi = i
180#                if dem_elevation_r[i,j] == NODATA_value:
181#                    nn += 1
182#
183#                if k == 0:
184#                    i1_0 = i
185#                    j1_0 = j
186#
187#                k += 1
188#
189#    index1 = j1_0
190#    index2 = thisj
191#
192#    # Dimension definitions
193#    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
194#    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
195#
196#    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
197#    nopoints = clippednopoints-nn
198#
199#    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
200#
201#    if verbose:
202#        log.critical('There are %d values in the elevation' % totalnopoints)
203#        log.critical('There are %d values in the clipped elevation'
204#                     % clippednopoints)
205#        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
206#
207#    outfile.createDimension('number_of_points', nopoints)
208#    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
209#
210#    # Variable definitions
211#    outfile.createVariable('points', netcdf_float, ('number_of_points',
212#                                                    'number_of_dimensions'))
213#    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
214#
215#    # Get handles to the variables
216#    points = outfile.variables['points']
217#    elevation = outfile.variables['elevation']
218#
219#    # Number of points
220#    N = points.shape[0]
221#
222#    lenv = index2-index1+1
223#
224#    # Store data
225#    global_index = 0
226#    # for i in range(nrows):
227#    for i in range(i1_0, thisi+1, 1):
228#        if verbose and i % ((nrows+10)/10) == 0:
229#            log.critical('Processing row %d of %d' % (i, nrows))
230#
231#        lower_index = global_index
232#
233#        v = dem_elevation_r[i,index1:index2+1]
234#        no_NODATA = num.sum(v == NODATA_value)
235#        if no_NODATA > 0:
236#            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
237#        else:
238#            newcols = lenv              # ncols_in_bounding_box
239#
240#        telev = num.zeros(newcols, num.float)
241#        tpoints = num.zeros((newcols, 2), num.float)
242#
243#        local_index = 0
244#
245#        y = (nrows-i-1)*cellsize + yllcorner
246#        #for j in range(ncols):
247#        for j in range(j1_0,index2+1,1):
248#            x = j*cellsize + xllcorner
249#            if easting_min <= x <= easting_max \
250#               and northing_min <= y <= northing_max \
251#               and dem_elevation_r[i,j] != NODATA_value:
252#
253#                #print [x-easting_min, y-northing_min]
254#                #print x , y
255#                #print easting_min, northing_min
256#                #print xllcorner, yllcorner
257#                #print cellsize
258#
259#                tpoints[local_index, :] = [x-easting_min, y-northing_min]
260#                telev[local_index] = dem_elevation_r[i, j]
261#                global_index += 1
262#                local_index += 1
263#
264#        upper_index = global_index
265#
266#        if upper_index == lower_index + newcols:
267#
268#            # Seems to be an error with the windows version of
269#            # Netcdf. The following gave errors
270#
271#            try:
272#                points[lower_index:upper_index, :] = tpoints
273#                elevation[lower_index:upper_index] = telev
274#            except:
275#                # so used the following if an error occurs
276#                for index in range(newcols):
277#                    points[index+lower_index, :] = tpoints[index,:]
278#                    elevation[index+lower_index] = telev[index]
279#
280#    assert global_index == nopoints, 'index not equal to number of points'
281
282
283    #========================================
284    # Do the preceeding with numpy
285    #========================================
286    y = num.arange(nrows,dtype=num.float)
287    y = yllcorner + (nrows-1)*cellsize - y*cellsize
288
289    x = num.arange(ncols,dtype=num.float)
290    x = xllcorner + x*cellsize
291
292    xx,yy = num.meshgrid(x,y)
293
294    xx = xx.flatten()
295    yy = yy.flatten()
296
297   
298    flag = num.logical_and(num.logical_and((xx <= easting_max),(xx >= easting_min)),
299                           num.logical_and((yy <= northing_max),(yy >= northing_min)))
300
301   
302    dem = dem_elevation[:].flatten()
303
304
305    id = num.where(flag)[0]
306
307    xx = xx[id]
308    yy = yy[id]
309    dem = dem[id]
310
311
312    clippednopoints = len(dem)
313    #print clippedpoints
314   
315    #print xx
316    #print yy
317    #print dem
318
319    data_flag = dem != NODATA_value
320
321    data_id = num.where(data_flag)
322
323    xx = xx[data_id]
324    yy = yy[data_id]
325    dem = dem[data_id]
326
327    nn = clippednopoints - len(dem)
328
329    nopoints = len(dem)
330
331
332    if verbose:
333        log.critical('There are %d values in the elevation' % totalnopoints)
334        log.critical('There are %d values in the clipped elevation'
335                     % clippednopoints)
336        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
337
338    outfile.createDimension('number_of_points', nopoints)
339    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
340
341    # Variable definitions
342    outfile.createVariable('points', netcdf_float, ('number_of_points',
343                                                    'number_of_dimensions'))
344    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
345
346    # Get handles to the variables
347    points = outfile.variables['points']
348    elevation = outfile.variables['elevation']
349
350    points[:,0] = xx - easting_min
351    points[:,1] = yy - northing_min
352    elevation[:] = dem
353
354
355    infile.close()
356    outfile.close()
357
Note: See TracBrowser for help on using the repository browser.