| 1 | # external modules |
|---|
| 2 | import numpy as num |
|---|
| 3 | |
|---|
| 4 | # ANUGA modules |
|---|
| 5 | import anuga.utilities.log as log |
|---|
| 6 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a, \ |
|---|
| 7 | netcdf_float |
|---|
| 8 | |
|---|
| 9 | from asc2dem import asc2dem |
|---|
| 10 | |
|---|
| 11 | |
|---|
| 12 | def 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 anuga.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 | |
|---|
| 55 | def _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 | |
|---|
