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

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

Further split up file conversion to modularise shallow_water module.

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