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

Last change on this file since 8466 was 8249, checked in by steve, 13 years ago

Got rid of those annoying double_scalar warnings in the windows code (just divide by zero warnings)
File size: 8.0 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 Scientific.IO.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 = infile.ncols[0]
89    nrows = infile.nrows[0]
90    xllcorner = infile.xllcorner[0]  # Easting of lower left corner
91    yllcorner = infile.yllcorner[0]  # Northing of lower left corner
92    cellsize = infile.cellsize[0]
93    NODATA_value = infile.NODATA_value[0]
94    dem_elevation = infile.variables['elevation']
95
96    zone = infile.zone[0]
97    false_easting = infile.false_easting[0]
98    false_northing = infile.false_northing[0]
99
100    # Text strings
101    projection = infile.projection
102    datum = infile.datum
103    units = infile.units
104
105    # Get output file
106    if name_out == None:
107        ptsname = root + '.pts'
108    else:
109        ptsname = name_out
110
111    if verbose: log.critical('Store to NetCDF file %s' % ptsname)
112
113    # NetCDF file definition
114    outfile = NetCDFFile(ptsname, netcdf_mode_w)
115
116    # Create new file
117    outfile.institution = 'Geoscience Australia'
118    outfile.description = 'NetCDF pts format for compact and portable ' \
119                          'storage of spatial point data'
120
121    # Assign default values
122    if easting_min is None: easting_min = xllcorner
123    if easting_max is None: easting_max = xllcorner + ncols*cellsize
124    if northing_min is None: northing_min = yllcorner
125    if northing_max is None: northing_max = yllcorner + nrows*cellsize
126
127    # Compute offsets to update georeferencing
128    easting_offset = xllcorner - easting_min
129    northing_offset = yllcorner - northing_min
130
131    # Georeferencing
132    outfile.zone = zone
133    outfile.xllcorner = easting_min # Easting of lower left corner
134    outfile.yllcorner = northing_min # Northing of lower left corner
135    outfile.false_easting = false_easting
136    outfile.false_northing = false_northing
137
138    outfile.projection = projection
139    outfile.datum = datum
140    outfile.units = units
141
142    # Grid info (FIXME: probably not going to be used, but heck)
143    outfile.ncols = ncols
144    outfile.nrows = nrows
145
146    dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols))
147    totalnopoints = nrows*ncols
148
149    # Calculating number of NODATA_values for each row in clipped region
150    # FIXME: use array operations to do faster
151    nn = 0
152    k = 0
153    i1_0 = 0
154    j1_0 = 0
155    thisj = 0
156    thisi = 0
157    for i in range(nrows):
158        y = (nrows-i-1)*cellsize + yllcorner
159        for j in range(ncols):
160            x = j*cellsize + xllcorner
161            if easting_min <= x <= easting_max \
162               and northing_min <= y <= northing_max:
163                thisj = j
164                thisi = i
165                if dem_elevation_r[i,j] == NODATA_value:
166                    nn += 1
167
168                if k == 0:
169                    i1_0 = i
170                    j1_0 = j
171
172                k += 1
173
174    index1 = j1_0
175    index2 = thisj
176
177    # Dimension definitions
178    nrows_in_bounding_box = int(round((northing_max-northing_min)/cellsize))
179    ncols_in_bounding_box = int(round((easting_max-easting_min)/cellsize))
180
181    clippednopoints = (thisi+1-i1_0)*(thisj+1-j1_0)
182    nopoints = clippednopoints-nn
183
184    clipped_dem_elev = dem_elevation_r[i1_0:thisi+1,j1_0:thisj+1]
185
186    if verbose:
187        log.critical('There are %d values in the elevation' % totalnopoints)
188        log.critical('There are %d values in the clipped elevation'
189                     % clippednopoints)
190        log.critical('There are %d NODATA_values in the clipped elevation' % nn)
191
192    outfile.createDimension('number_of_points', nopoints)
193    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
194
195    # Variable definitions
196    outfile.createVariable('points', netcdf_float, ('number_of_points',
197                                                    'number_of_dimensions'))
198    outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
199
200    # Get handles to the variables
201    points = outfile.variables['points']
202    elevation = outfile.variables['elevation']
203
204    lenv = index2-index1+1
205
206    # Store data
207    global_index = 0
208    # for i in range(nrows):
209    for i in range(i1_0, thisi+1, 1):
210        if verbose and i % ((nrows+10)/10) == 0:
211            log.critical('Processing row %d of %d' % (i, nrows))
212
213        lower_index = global_index
214
215        v = dem_elevation_r[i,index1:index2+1]
216        no_NODATA = num.sum(v == NODATA_value)
217        if no_NODATA > 0:
218            newcols = lenv - no_NODATA  # ncols_in_bounding_box - no_NODATA
219        else:
220            newcols = lenv              # ncols_in_bounding_box
221
222        telev = num.zeros(newcols, num.float)
223        tpoints = num.zeros((newcols, 2), num.float)
224
225        local_index = 0
226
227        y = (nrows-i-1)*cellsize + yllcorner
228        #for j in range(ncols):
229        for j in range(j1_0,index2+1,1):
230            x = j*cellsize + xllcorner
231            if easting_min <= x <= easting_max \
232               and northing_min <= y <= northing_max \
233               and dem_elevation_r[i,j] != NODATA_value:
234                tpoints[local_index, :] = [x-easting_min, y-northing_min]
235                telev[local_index] = dem_elevation_r[i, j]
236                global_index += 1
237                local_index += 1
238
239        upper_index = global_index
240
241        if upper_index == lower_index + newcols:
242            # Seems to be an error with the windows version of
243            # Netcdf. The following gave errors
244            #points[lower_index:upper_index, :] = tpoints
245            #elevation[lower_index:upper_index] = telev
246            # so used the following
247            for index in range(newcols):
248                points[index+lower_index, :] = tpoints[index,:]
249                elevation[index+lower_index] = telev[index]
250
251    assert global_index == nopoints, 'index not equal to number of points'
252
253    infile.close()
254    outfile.close()
255
Note: See TracBrowser for help on using the repository browser.