| 1 | """ Load a DEM file, decimate it, and resave it. |
|---|
| 2 | """ |
|---|
| 3 | |
|---|
| 4 | import numpy as num |
|---|
| 5 | |
|---|
| 6 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_float |
|---|
| 7 | |
|---|
| 8 | def dem2dem(basename_in, stencil, cellsize_new, basename_out=None, |
|---|
| 9 | verbose=False): |
|---|
| 10 | """Read Digitial Elevation model from the following NetCDF format (.dem) |
|---|
| 11 | |
|---|
| 12 | Example: |
|---|
| 13 | |
|---|
| 14 | ncols 3121 |
|---|
| 15 | nrows 1800 |
|---|
| 16 | xllcorner 722000 |
|---|
| 17 | yllcorner 5893000 |
|---|
| 18 | cellsize 25 |
|---|
| 19 | NODATA_value -9999 |
|---|
| 20 | 138.3698 137.4194 136.5062 135.5558 .......... |
|---|
| 21 | |
|---|
| 22 | Decimate data to cellsize_new using stencil and write to NetCDF dem format. |
|---|
| 23 | """ |
|---|
| 24 | |
|---|
| 25 | import os |
|---|
| 26 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 27 | |
|---|
| 28 | root = basename_in |
|---|
| 29 | inname = root + '.dem' |
|---|
| 30 | |
|---|
| 31 | #Open existing netcdf file to read |
|---|
| 32 | infile = NetCDFFile(inname, netcdf_mode_r) |
|---|
| 33 | |
|---|
| 34 | if verbose: log.critical('Reading DEM from %s' % inname) |
|---|
| 35 | |
|---|
| 36 | # Read metadata (convert from numpy.int32 to int where appropriate) |
|---|
| 37 | ncols = int(infile.ncols[0]) |
|---|
| 38 | nrows = int(infile.nrows[0]) |
|---|
| 39 | xllcorner = infile.xllcorner[0] |
|---|
| 40 | yllcorner = infile.yllcorner[0] |
|---|
| 41 | cellsize = int(infile.cellsize[0]) |
|---|
| 42 | NODATA_value = int(infile.NODATA_value[0]) |
|---|
| 43 | zone = int(infile.zone[0]) |
|---|
| 44 | false_easting = infile.false_easting[0] |
|---|
| 45 | false_northing = infile.false_northing[0] |
|---|
| 46 | projection = infile.projection |
|---|
| 47 | datum = infile.datum |
|---|
| 48 | units = infile.units |
|---|
| 49 | |
|---|
| 50 | dem_elevation = infile.variables['elevation'] |
|---|
| 51 | |
|---|
| 52 | #Get output file name |
|---|
| 53 | if basename_out == None: |
|---|
| 54 | outname = root + '_' + repr(cellsize_new) + '.dem' |
|---|
| 55 | else: |
|---|
| 56 | outname = basename_out + '.dem' |
|---|
| 57 | |
|---|
| 58 | if verbose: log.critical('Write decimated NetCDF file to %s' % outname) |
|---|
| 59 | |
|---|
| 60 | #Determine some dimensions for decimated grid |
|---|
| 61 | (nrows_stencil, ncols_stencil) = stencil.shape |
|---|
| 62 | x_offset = ncols_stencil / 2 |
|---|
| 63 | y_offset = nrows_stencil / 2 |
|---|
| 64 | cellsize_ratio = int(cellsize_new / cellsize) |
|---|
| 65 | ncols_new = 1 + (ncols - ncols_stencil) / cellsize_ratio |
|---|
| 66 | nrows_new = 1 + (nrows - nrows_stencil) / cellsize_ratio |
|---|
| 67 | |
|---|
| 68 | #print type(ncols_new), ncols_new |
|---|
| 69 | |
|---|
| 70 | #Open netcdf file for output |
|---|
| 71 | outfile = NetCDFFile(outname, netcdf_mode_w) |
|---|
| 72 | |
|---|
| 73 | #Create new file |
|---|
| 74 | outfile.institution = 'Geoscience Australia' |
|---|
| 75 | outfile.description = 'NetCDF DEM format for compact and portable ' \ |
|---|
| 76 | 'storage of spatial point data' |
|---|
| 77 | |
|---|
| 78 | #Georeferencing |
|---|
| 79 | outfile.zone = zone |
|---|
| 80 | outfile.projection = projection |
|---|
| 81 | outfile.datum = datum |
|---|
| 82 | outfile.units = units |
|---|
| 83 | |
|---|
| 84 | outfile.cellsize = cellsize_new |
|---|
| 85 | outfile.NODATA_value = NODATA_value |
|---|
| 86 | outfile.false_easting = false_easting |
|---|
| 87 | outfile.false_northing = false_northing |
|---|
| 88 | |
|---|
| 89 | outfile.xllcorner = xllcorner + (x_offset * cellsize) |
|---|
| 90 | outfile.yllcorner = yllcorner + (y_offset * cellsize) |
|---|
| 91 | outfile.ncols = ncols_new |
|---|
| 92 | outfile.nrows = nrows_new |
|---|
| 93 | |
|---|
| 94 | # dimension definition |
|---|
| 95 | #print nrows_new, ncols_new, nrows_new*ncols_new |
|---|
| 96 | #print type(nrows_new), type(ncols_new), type(nrows_new*ncols_new) |
|---|
| 97 | outfile.createDimension('number_of_points', nrows_new*ncols_new) |
|---|
| 98 | |
|---|
| 99 | # variable definition |
|---|
| 100 | outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) |
|---|
| 101 | |
|---|
| 102 | # Get handle to the variable |
|---|
| 103 | elevation = outfile.variables['elevation'] |
|---|
| 104 | |
|---|
| 105 | dem_elevation_r = num.reshape(dem_elevation, (nrows, ncols)) |
|---|
| 106 | |
|---|
| 107 | #Store data |
|---|
| 108 | global_index = 0 |
|---|
| 109 | for i in range(nrows_new): |
|---|
| 110 | if verbose: log.critical('Processing row %d of %d' % (i, nrows_new)) |
|---|
| 111 | |
|---|
| 112 | lower_index = global_index |
|---|
| 113 | telev = num.zeros(ncols_new, num.float) |
|---|
| 114 | local_index = 0 |
|---|
| 115 | trow = i * cellsize_ratio |
|---|
| 116 | |
|---|
| 117 | for j in range(ncols_new): |
|---|
| 118 | tcol = j * cellsize_ratio |
|---|
| 119 | tmp = dem_elevation_r[trow:trow+nrows_stencil, |
|---|
| 120 | tcol:tcol+ncols_stencil] |
|---|
| 121 | |
|---|
| 122 | #if dem contains 1 or more NODATA_values set value in |
|---|
| 123 | #decimated dem to NODATA_value, else compute decimated |
|---|
| 124 | #value using stencil |
|---|
| 125 | if num.sum(num.sum(num.equal(tmp, NODATA_value))) > 0: |
|---|
| 126 | telev[local_index] = NODATA_value |
|---|
| 127 | else: |
|---|
| 128 | telev[local_index] = num.sum(num.sum(tmp * stencil)) |
|---|
| 129 | |
|---|
| 130 | global_index += 1 |
|---|
| 131 | local_index += 1 |
|---|
| 132 | |
|---|
| 133 | upper_index = global_index |
|---|
| 134 | |
|---|
| 135 | elevation[lower_index:upper_index] = telev |
|---|
| 136 | |
|---|
| 137 | assert global_index == nrows_new*ncols_new, \ |
|---|
| 138 | 'index not equal to number of points' |
|---|
| 139 | |
|---|
| 140 | infile.close() |
|---|
| 141 | outfile.close() |
|---|
| 142 | |
|---|
