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

Last change on this file since 7804 was 7804, checked in by hudson, 12 years ago

Fixed up failing tests, updated user guide with new API (first few chapters only).

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