[2344] | 1 | |
---|
| 2 | |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | ############################################### |
---|
| 9 | #OBSOLETE STUFF |
---|
| 10 | #Native checkpoint format. |
---|
| 11 | #Information needed to recreate a state is preserved |
---|
| 12 | #FIXME: Rethink and maybe use netcdf format |
---|
| 13 | def cpt_variable_writer(filename, t, v0, v1, v2): |
---|
| 14 | """Store all conserved quantities to file |
---|
| 15 | """ |
---|
| 16 | |
---|
| 17 | M, N = v0.shape |
---|
| 18 | |
---|
| 19 | FN = create_filename(filename, 'cpt', M, t) |
---|
| 20 | #print 'Writing to %s' %FN |
---|
| 21 | |
---|
| 22 | fid = open(FN, 'w') |
---|
| 23 | for i in range(M): |
---|
| 24 | for j in range(N): |
---|
| 25 | fid.write('%.16e ' %v0[i,j]) |
---|
| 26 | for j in range(N): |
---|
| 27 | fid.write('%.16e ' %v1[i,j]) |
---|
| 28 | for j in range(N): |
---|
| 29 | fid.write('%.16e ' %v2[i,j]) |
---|
| 30 | |
---|
| 31 | fid.write('\n') |
---|
| 32 | fid.close() |
---|
| 33 | |
---|
| 34 | |
---|
| 35 | def cpt_variable_reader(filename, t, v0, v1, v2): |
---|
| 36 | """Store all conserved quantities to file |
---|
| 37 | """ |
---|
| 38 | |
---|
| 39 | M, N = v0.shape |
---|
| 40 | |
---|
| 41 | FN = create_filename(filename, 'cpt', M, t) |
---|
| 42 | #print 'Reading from %s' %FN |
---|
| 43 | |
---|
| 44 | fid = open(FN) |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | for i in range(M): |
---|
| 48 | values = fid.readline().split() #Get one line |
---|
| 49 | |
---|
| 50 | for j in range(N): |
---|
| 51 | v0[i,j] = float(values[j]) |
---|
| 52 | v1[i,j] = float(values[3+j]) |
---|
| 53 | v2[i,j] = float(values[6+j]) |
---|
| 54 | |
---|
| 55 | fid.close() |
---|
| 56 | |
---|
| 57 | def cpt_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
| 58 | """Writes x,y,z,z,z coordinates of triangles constituting the bed |
---|
| 59 | elevation. |
---|
| 60 | FIXME: Not in use pt |
---|
| 61 | """ |
---|
| 62 | |
---|
| 63 | M, N = v0.shape |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | print X0 |
---|
| 67 | import sys; sys.exit() |
---|
| 68 | FN = create_filename(filename, 'cpt', M) |
---|
| 69 | print 'Writing to %s' %FN |
---|
| 70 | |
---|
| 71 | fid = open(FN, 'w') |
---|
| 72 | for i in range(M): |
---|
| 73 | for j in range(2): |
---|
| 74 | fid.write('%.16e ' %X0[i,j]) #x, y |
---|
| 75 | for j in range(N): |
---|
| 76 | fid.write('%.16e ' %v0[i,j]) #z,z,z, |
---|
| 77 | |
---|
| 78 | for j in range(2): |
---|
| 79 | fid.write('%.16e ' %X1[i,j]) #x, y |
---|
| 80 | for j in range(N): |
---|
| 81 | fid.write('%.16e ' %v1[i,j]) |
---|
| 82 | |
---|
| 83 | for j in range(2): |
---|
| 84 | fid.write('%.16e ' %X2[i,j]) #x, y |
---|
| 85 | for j in range(N): |
---|
| 86 | fid.write('%.16e ' %v2[i,j]) |
---|
| 87 | |
---|
| 88 | fid.write('\n') |
---|
| 89 | fid.close() |
---|
| 90 | |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | #Function for storing out to e.g. visualisation |
---|
| 94 | #FIXME: Do we want this? |
---|
| 95 | #FIXME: Not done yet for this version |
---|
| 96 | def dat_constant_writer(filename, X0, X1, X2, v0, v1, v2): |
---|
| 97 | """Writes x,y,z coordinates of triangles constituting the bed elevation. |
---|
| 98 | """ |
---|
| 99 | |
---|
| 100 | M, N = v0.shape |
---|
| 101 | |
---|
| 102 | FN = create_filename(filename, 'dat', M) |
---|
| 103 | #print 'Writing to %s' %FN |
---|
| 104 | |
---|
| 105 | fid = open(FN, 'w') |
---|
| 106 | for i in range(M): |
---|
| 107 | for j in range(2): |
---|
| 108 | fid.write('%f ' %X0[i,j]) #x, y |
---|
| 109 | fid.write('%f ' %v0[i,0]) #z |
---|
| 110 | |
---|
| 111 | for j in range(2): |
---|
| 112 | fid.write('%f ' %X1[i,j]) #x, y |
---|
| 113 | fid.write('%f ' %v1[i,0]) #z |
---|
| 114 | |
---|
| 115 | for j in range(2): |
---|
| 116 | fid.write('%f ' %X2[i,j]) #x, y |
---|
| 117 | fid.write('%f ' %v2[i,0]) #z |
---|
| 118 | |
---|
| 119 | fid.write('\n') |
---|
| 120 | fid.close() |
---|
| 121 | |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | def dat_variable_writer(filename, t, v0, v1, v2): |
---|
| 125 | """Store water height to file |
---|
| 126 | """ |
---|
| 127 | |
---|
| 128 | M, N = v0.shape |
---|
| 129 | |
---|
| 130 | FN = create_filename(filename, 'dat', M, t) |
---|
| 131 | #print 'Writing to %s' %FN |
---|
| 132 | |
---|
| 133 | fid = open(FN, 'w') |
---|
| 134 | for i in range(M): |
---|
| 135 | fid.write('%.4f ' %v0[i,0]) |
---|
| 136 | fid.write('%.4f ' %v1[i,0]) |
---|
| 137 | fid.write('%.4f ' %v2[i,0]) |
---|
| 138 | |
---|
| 139 | fid.write('\n') |
---|
| 140 | fid.close() |
---|
| 141 | |
---|
| 142 | |
---|
| 143 | def read_sww(filename): |
---|
| 144 | """Read sww Net CDF file containing Shallow Water Wave simulation |
---|
| 145 | |
---|
| 146 | The integer array volumes is of shape Nx3 where N is the number of |
---|
| 147 | triangles in the mesh. |
---|
| 148 | |
---|
| 149 | Each entry in volumes is an index into the x,y arrays (the location). |
---|
| 150 | |
---|
| 151 | Quantities stage, elevation, xmomentum and ymomentum are all in arrays of dimensions |
---|
| 152 | number_of_timesteps, number_of_points. |
---|
| 153 | |
---|
| 154 | The momentum is not always stored. |
---|
| 155 | |
---|
| 156 | """ |
---|
| 157 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 158 | print 'Reading from ', filename |
---|
| 159 | fid = NetCDFFile(filename, 'r') #Open existing file for read |
---|
| 160 | #latitude, longitude |
---|
| 161 | # Get the variables as Numeric arrays |
---|
| 162 | x = fid.variables['x'] #x-coordinates of vertices |
---|
| 163 | y = fid.variables['y'] #y-coordinates of vertices |
---|
| 164 | z = fid.variables['elevation'] #Elevation |
---|
| 165 | time = fid.variables['time'] #Timesteps |
---|
| 166 | stage = fid.variables['stage'] #Water level |
---|
| 167 | #xmomentum = fid.variables['xmomentum'] #Momentum in the x-direction |
---|
| 168 | #ymomentum = fid.variables['ymomentum'] #Momentum in the y-direction |
---|
| 169 | |
---|
| 170 | volumes = fid.variables['volumes'] #Connectivity |
---|
| 171 | |
---|
| 172 | #FIXME (Ole): What is this? |
---|
| 173 | # Why isn't anything returned? |
---|
| 174 | # Where's the unit test? |
---|
| 175 | |
---|
| 176 | |
---|
| 177 | |
---|
| 178 | def sww2asc_obsolete(basename_in, basename_out = None, |
---|
| 179 | quantity = None, |
---|
| 180 | timestep = None, |
---|
| 181 | reduction = None, |
---|
| 182 | cellsize = 10, |
---|
| 183 | verbose = False, |
---|
| 184 | origin = None): |
---|
| 185 | """Read SWW file and convert to Digitial Elevation model format (.asc) |
---|
| 186 | |
---|
| 187 | Example: |
---|
| 188 | |
---|
| 189 | ncols 3121 |
---|
| 190 | nrows 1800 |
---|
| 191 | xllcorner 722000 |
---|
| 192 | yllcorner 5893000 |
---|
| 193 | cellsize 25 |
---|
| 194 | NODATA_value -9999 |
---|
| 195 | 138.3698 137.4194 136.5062 135.5558 .......... |
---|
| 196 | |
---|
| 197 | Also write accompanying file with same basename_in but extension .prj |
---|
| 198 | used to fix the UTM zone, datum, false northings and eastings. |
---|
| 199 | |
---|
| 200 | The prj format is assumed to be as |
---|
| 201 | |
---|
| 202 | Projection UTM |
---|
| 203 | Zone 56 |
---|
| 204 | Datum WGS84 |
---|
| 205 | Zunits NO |
---|
| 206 | Units METERS |
---|
| 207 | Spheroid WGS84 |
---|
| 208 | Xshift 0.0000000000 |
---|
| 209 | Yshift 10000000.0000000000 |
---|
| 210 | Parameters |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | if quantity is given, out values from quantity otherwise default to |
---|
| 214 | elevation |
---|
| 215 | |
---|
| 216 | if timestep (an index) is given, output quantity at that timestep |
---|
| 217 | |
---|
| 218 | if reduction is given use that to reduce quantity over all timesteps. |
---|
| 219 | |
---|
| 220 | """ |
---|
| 221 | from Numeric import array, Float, concatenate, NewAxis, zeros,\ |
---|
| 222 | sometrue |
---|
[3514] | 223 | from anuga.utilities.polygon import inside_polygon |
---|
[2344] | 224 | |
---|
| 225 | #FIXME: Should be variable |
---|
| 226 | datum = 'WGS84' |
---|
| 227 | false_easting = 500000 |
---|
| 228 | false_northing = 10000000 |
---|
| 229 | |
---|
| 230 | if quantity is None: |
---|
| 231 | quantity = 'elevation' |
---|
| 232 | |
---|
| 233 | if reduction is None: |
---|
| 234 | reduction = max |
---|
| 235 | |
---|
| 236 | if basename_out is None: |
---|
| 237 | basename_out = basename_in + '_%s' %quantity |
---|
| 238 | |
---|
| 239 | swwfile = basename_in + '.sww' |
---|
| 240 | ascfile = basename_out + '.asc' |
---|
| 241 | prjfile = basename_out + '.prj' |
---|
| 242 | |
---|
| 243 | |
---|
| 244 | if verbose: print 'Reading from %s' %swwfile |
---|
| 245 | #Read sww file |
---|
| 246 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 247 | fid = NetCDFFile(swwfile) |
---|
| 248 | |
---|
| 249 | #Get extent and reference |
---|
| 250 | x = fid.variables['x'][:] |
---|
| 251 | y = fid.variables['y'][:] |
---|
| 252 | volumes = fid.variables['volumes'][:] |
---|
| 253 | |
---|
| 254 | ymin = min(y); ymax = max(y) |
---|
| 255 | xmin = min(x); xmax = max(x) |
---|
| 256 | |
---|
| 257 | number_of_timesteps = fid.dimensions['number_of_timesteps'] |
---|
| 258 | number_of_points = fid.dimensions['number_of_points'] |
---|
| 259 | if origin is None: |
---|
| 260 | |
---|
| 261 | #Get geo_reference |
---|
| 262 | #sww files don't have to have a geo_ref |
---|
| 263 | try: |
---|
| 264 | geo_reference = Geo_reference(NetCDFObject=fid) |
---|
| 265 | except AttributeError, e: |
---|
| 266 | geo_reference = Geo_reference() #Default georef object |
---|
| 267 | |
---|
| 268 | xllcorner = geo_reference.get_xllcorner() |
---|
| 269 | yllcorner = geo_reference.get_yllcorner() |
---|
| 270 | zone = geo_reference.get_zone() |
---|
| 271 | else: |
---|
| 272 | zone = origin[0] |
---|
| 273 | xllcorner = origin[1] |
---|
| 274 | yllcorner = origin[2] |
---|
| 275 | |
---|
| 276 | |
---|
| 277 | #Get quantity and reduce if applicable |
---|
| 278 | if verbose: print 'Reading quantity %s' %quantity |
---|
| 279 | |
---|
| 280 | if quantity.lower() == 'depth': |
---|
| 281 | q = fid.variables['stage'][:] - fid.variables['elevation'][:] |
---|
| 282 | else: |
---|
| 283 | q = fid.variables[quantity][:] |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | if len(q.shape) == 2: |
---|
| 287 | if verbose: print 'Reducing quantity %s' %quantity |
---|
| 288 | q_reduced = zeros( number_of_points, Float ) |
---|
| 289 | |
---|
| 290 | for k in range(number_of_points): |
---|
| 291 | q_reduced[k] = reduction( q[:,k] ) |
---|
| 292 | |
---|
| 293 | q = q_reduced |
---|
| 294 | |
---|
| 295 | #Now q has dimension: number_of_points |
---|
| 296 | |
---|
| 297 | #Create grid and update xll/yll corner and x,y |
---|
| 298 | if verbose: print 'Creating grid' |
---|
| 299 | ncols = int((xmax-xmin)/cellsize)+1 |
---|
| 300 | nrows = int((ymax-ymin)/cellsize)+1 |
---|
| 301 | |
---|
| 302 | newxllcorner = xmin+xllcorner |
---|
| 303 | newyllcorner = ymin+yllcorner |
---|
| 304 | |
---|
| 305 | x = x+xllcorner-newxllcorner |
---|
| 306 | y = y+yllcorner-newyllcorner |
---|
| 307 | |
---|
| 308 | vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1) |
---|
| 309 | assert len(vertex_points.shape) == 2 |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | from Numeric import zeros, Float |
---|
| 313 | grid_points = zeros ( (ncols*nrows, 2), Float ) |
---|
| 314 | |
---|
| 315 | |
---|
| 316 | for i in xrange(nrows): |
---|
| 317 | yg = i*cellsize |
---|
| 318 | for j in xrange(ncols): |
---|
| 319 | xg = j*cellsize |
---|
| 320 | k = i*ncols + j |
---|
| 321 | |
---|
| 322 | grid_points[k,0] = xg |
---|
| 323 | grid_points[k,1] = yg |
---|
| 324 | |
---|
| 325 | #Interpolate |
---|
| 326 | from least_squares import Interpolation |
---|
| 327 | |
---|
| 328 | |
---|
| 329 | #FIXME: This should be done with precrop = True, otherwise it'll |
---|
| 330 | #take forever. With expand_search set to False, some grid points might |
---|
| 331 | #miss out.... |
---|
| 332 | |
---|
| 333 | interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0, |
---|
| 334 | precrop = False, expand_search = True, |
---|
| 335 | verbose = verbose) |
---|
| 336 | |
---|
| 337 | #Interpolate using quantity values |
---|
| 338 | if verbose: print 'Interpolating' |
---|
| 339 | grid_values = interp.interpolate(q).flat |
---|
| 340 | |
---|
| 341 | #Write |
---|
| 342 | #Write prj file |
---|
| 343 | if verbose: print 'Writing %s' %prjfile |
---|
| 344 | prjid = open(prjfile, 'w') |
---|
| 345 | prjid.write('Projection %s\n' %'UTM') |
---|
| 346 | prjid.write('Zone %d\n' %zone) |
---|
| 347 | prjid.write('Datum %s\n' %datum) |
---|
| 348 | prjid.write('Zunits NO\n') |
---|
| 349 | prjid.write('Units METERS\n') |
---|
| 350 | prjid.write('Spheroid %s\n' %datum) |
---|
| 351 | prjid.write('Xshift %d\n' %false_easting) |
---|
| 352 | prjid.write('Yshift %d\n' %false_northing) |
---|
| 353 | prjid.write('Parameters\n') |
---|
| 354 | prjid.close() |
---|
| 355 | |
---|
| 356 | |
---|
| 357 | |
---|
| 358 | if verbose: print 'Writing %s' %ascfile |
---|
| 359 | NODATA_value = -9999 |
---|
| 360 | |
---|
| 361 | ascid = open(ascfile, 'w') |
---|
| 362 | |
---|
| 363 | ascid.write('ncols %d\n' %ncols) |
---|
| 364 | ascid.write('nrows %d\n' %nrows) |
---|
| 365 | ascid.write('xllcorner %d\n' %newxllcorner) |
---|
| 366 | ascid.write('yllcorner %d\n' %newyllcorner) |
---|
| 367 | ascid.write('cellsize %f\n' %cellsize) |
---|
| 368 | ascid.write('NODATA_value %d\n' %NODATA_value) |
---|
| 369 | |
---|
| 370 | |
---|
| 371 | #Get bounding polygon from mesh |
---|
| 372 | P = interp.mesh.get_boundary_polygon() |
---|
| 373 | inside_indices = inside_polygon(grid_points, P) |
---|
| 374 | |
---|
| 375 | for i in range(nrows): |
---|
| 376 | if verbose and i%((nrows+10)/10)==0: |
---|
| 377 | print 'Doing row %d of %d' %(i, nrows) |
---|
| 378 | |
---|
| 379 | for j in range(ncols): |
---|
| 380 | index = (nrows-i-1)*ncols+j |
---|
| 381 | |
---|
| 382 | if sometrue(inside_indices == index): |
---|
| 383 | ascid.write('%f ' %grid_values[index]) |
---|
| 384 | else: |
---|
| 385 | ascid.write('%d ' %NODATA_value) |
---|
| 386 | |
---|
| 387 | ascid.write('\n') |
---|
| 388 | |
---|
| 389 | #Close |
---|
| 390 | ascid.close() |
---|
| 391 | fid.close() |
---|
| 392 | |
---|
[3761] | 393 | |
---|
| 394 | |
---|
| 395 | |
---|
| 396 | def NOT_test_dem2pts(self): |
---|
| 397 | """Test conversion from dem in ascii format to native NetCDF xya format |
---|
| 398 | """ |
---|
| 399 | |
---|
| 400 | import time, os |
---|
| 401 | from Numeric import array, zeros, allclose, Float, concatenate |
---|
| 402 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 403 | |
---|
| 404 | #Write test asc file |
---|
| 405 | root = 'demtest' |
---|
| 406 | |
---|
| 407 | filename = root+'.asc' |
---|
| 408 | fid = open(filename, 'w') |
---|
| 409 | fid.write("""ncols 5 |
---|
| 410 | nrows 6 |
---|
| 411 | xllcorner 2000.5 |
---|
| 412 | yllcorner 3000.5 |
---|
| 413 | cellsize 25 |
---|
| 414 | NODATA_value -9999 |
---|
| 415 | """) |
---|
| 416 | #Create linear function |
---|
| 417 | |
---|
| 418 | ref_points = [] |
---|
| 419 | ref_elevation = [] |
---|
| 420 | for i in range(6): |
---|
| 421 | y = (6-i)*25.0 |
---|
| 422 | for j in range(5): |
---|
| 423 | x = j*25.0 |
---|
| 424 | z = x+2*y |
---|
| 425 | |
---|
| 426 | ref_points.append( [x,y] ) |
---|
| 427 | ref_elevation.append(z) |
---|
| 428 | fid.write('%f ' %z) |
---|
| 429 | fid.write('\n') |
---|
| 430 | |
---|
| 431 | fid.close() |
---|
| 432 | |
---|
| 433 | #Write prj file with metadata |
---|
| 434 | metafilename = root+'.prj' |
---|
| 435 | fid = open(metafilename, 'w') |
---|
| 436 | |
---|
| 437 | |
---|
| 438 | fid.write("""Projection UTM |
---|
| 439 | Zone 56 |
---|
| 440 | Datum WGS84 |
---|
| 441 | Zunits NO |
---|
| 442 | Units METERS |
---|
| 443 | Spheroid WGS84 |
---|
| 444 | Xshift 0.0000000000 |
---|
| 445 | Yshift 10000000.0000000000 |
---|
| 446 | Parameters |
---|
| 447 | """) |
---|
| 448 | fid.close() |
---|
| 449 | |
---|
| 450 | #Convert to NetCDF pts |
---|
| 451 | convert_dem_from_ascii2netcdf(root) |
---|
| 452 | dem2pts(root) |
---|
| 453 | |
---|
| 454 | #Check contents |
---|
| 455 | #Get NetCDF |
---|
| 456 | fid = NetCDFFile(root+'.pts', 'r') |
---|
| 457 | |
---|
| 458 | # Get the variables |
---|
| 459 | #print fid.variables.keys() |
---|
| 460 | points = fid.variables['points'] |
---|
| 461 | elevation = fid.variables['elevation'] |
---|
| 462 | |
---|
| 463 | #Check values |
---|
| 464 | |
---|
| 465 | #print points[:] |
---|
| 466 | #print ref_points |
---|
| 467 | assert allclose(points, ref_points) |
---|
| 468 | |
---|
| 469 | #print attributes[:] |
---|
| 470 | #print ref_elevation |
---|
| 471 | assert allclose(elevation, ref_elevation) |
---|
| 472 | |
---|
| 473 | #Cleanup |
---|
| 474 | fid.close() |
---|
| 475 | |
---|
| 476 | |
---|
| 477 | os.remove(root + '.pts') |
---|
| 478 | os.remove(root + '.dem') |
---|
| 479 | os.remove(root + '.asc') |
---|
| 480 | os.remove(root + '.prj') |
---|
| 481 | |
---|
| 482 | |
---|
| 483 | |
---|
| 484 | def NOT_test_dem2pts_bounding_box(self): |
---|
| 485 | """Test conversion from dem in ascii format to native NetCDF xya format |
---|
| 486 | """ |
---|
| 487 | |
---|
| 488 | import time, os |
---|
| 489 | from Numeric import array, zeros, allclose, Float, concatenate |
---|
| 490 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 491 | |
---|
| 492 | #Write test asc file |
---|
| 493 | root = 'demtest' |
---|
| 494 | |
---|
| 495 | filename = root+'.asc' |
---|
| 496 | fid = open(filename, 'w') |
---|
| 497 | fid.write("""ncols 5 |
---|
| 498 | nrows 6 |
---|
| 499 | xllcorner 2000.5 |
---|
| 500 | yllcorner 3000.5 |
---|
| 501 | cellsize 25 |
---|
| 502 | NODATA_value -9999 |
---|
| 503 | """) |
---|
| 504 | #Create linear function |
---|
| 505 | |
---|
| 506 | ref_points = [] |
---|
| 507 | ref_elevation = [] |
---|
| 508 | for i in range(6): |
---|
| 509 | y = (6-i)*25.0 |
---|
| 510 | for j in range(5): |
---|
| 511 | x = j*25.0 |
---|
| 512 | z = x+2*y |
---|
| 513 | |
---|
| 514 | ref_points.append( [x,y] ) |
---|
| 515 | ref_elevation.append(z) |
---|
| 516 | fid.write('%f ' %z) |
---|
| 517 | fid.write('\n') |
---|
| 518 | |
---|
| 519 | fid.close() |
---|
| 520 | |
---|
| 521 | #Write prj file with metadata |
---|
| 522 | metafilename = root+'.prj' |
---|
| 523 | fid = open(metafilename, 'w') |
---|
| 524 | |
---|
| 525 | |
---|
| 526 | fid.write("""Projection UTM |
---|
| 527 | Zone 56 |
---|
| 528 | Datum WGS84 |
---|
| 529 | Zunits NO |
---|
| 530 | Units METERS |
---|
| 531 | Spheroid WGS84 |
---|
| 532 | Xshift 0.0000000000 |
---|
| 533 | Yshift 10000000.0000000000 |
---|
| 534 | Parameters |
---|
| 535 | """) |
---|
| 536 | fid.close() |
---|
| 537 | |
---|
| 538 | #Convert to NetCDF pts |
---|
| 539 | convert_dem_from_ascii2netcdf(root) |
---|
| 540 | dem2pts(root, easting_min=2010.0, easting_max=2110.0, |
---|
| 541 | northing_min=3035.0, northing_max=3125.5) |
---|
| 542 | |
---|
| 543 | #Check contents |
---|
| 544 | #Get NetCDF |
---|
| 545 | fid = NetCDFFile(root+'.pts', 'r') |
---|
| 546 | |
---|
| 547 | # Get the variables |
---|
| 548 | #print fid.variables.keys() |
---|
| 549 | points = fid.variables['points'] |
---|
| 550 | elevation = fid.variables['elevation'] |
---|
| 551 | |
---|
| 552 | #Check values |
---|
| 553 | assert fid.xllcorner[0] == 2010.0 |
---|
| 554 | assert fid.yllcorner[0] == 3035.0 |
---|
| 555 | |
---|
| 556 | #create new reference points |
---|
| 557 | ref_points = [] |
---|
| 558 | ref_elevation = [] |
---|
| 559 | for i in range(4): |
---|
| 560 | y = (4-i)*25.0 + 25.0 |
---|
| 561 | y_new = y + 3000.5 - 3035.0 |
---|
| 562 | for j in range(4): |
---|
| 563 | x = j*25.0 + 25.0 |
---|
| 564 | x_new = x + 2000.5 - 2010.0 |
---|
| 565 | z = x+2*y |
---|
| 566 | |
---|
| 567 | ref_points.append( [x_new,y_new] ) |
---|
| 568 | ref_elevation.append(z) |
---|
| 569 | |
---|
| 570 | #print points[:] |
---|
| 571 | #print ref_points |
---|
| 572 | assert allclose(points, ref_points) |
---|
| 573 | |
---|
| 574 | #print attributes[:] |
---|
| 575 | #print ref_elevation |
---|
| 576 | assert allclose(elevation, ref_elevation) |
---|
| 577 | |
---|
| 578 | #Cleanup |
---|
| 579 | fid.close() |
---|
| 580 | |
---|
| 581 | |
---|
| 582 | os.remove(root + '.pts') |
---|
| 583 | os.remove(root + '.dem') |
---|
| 584 | os.remove(root + '.asc') |
---|
| 585 | os.remove(root + '.prj') |
---|
| 586 | |
---|
| 587 | |
---|
| 588 | |
---|
| 589 | def NOT_test_dem2pts_remove_Nullvalues(self): |
---|
| 590 | """Test conversion from dem in ascii format to native NetCDF xya format |
---|
| 591 | """ |
---|
| 592 | |
---|
| 593 | import time, os |
---|
| 594 | from Numeric import array, zeros, allclose, Float, concatenate |
---|
| 595 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 596 | |
---|
| 597 | #Write test asc file |
---|
| 598 | root = 'demtest' |
---|
| 599 | |
---|
| 600 | filename = root+'.asc' |
---|
| 601 | fid = open(filename, 'w') |
---|
| 602 | fid.write("""ncols 5 |
---|
| 603 | nrows 6 |
---|
| 604 | xllcorner 2000.5 |
---|
| 605 | yllcorner 3000.5 |
---|
| 606 | cellsize 25 |
---|
| 607 | NODATA_value -9999 |
---|
| 608 | """) |
---|
| 609 | #Create linear function |
---|
| 610 | # ref_ will write all the values |
---|
| 611 | # new_ref_ will write the values except for NODATA_values |
---|
| 612 | ref_points = [] |
---|
| 613 | ref_elevation = [] |
---|
| 614 | new_ref_pts = [] |
---|
| 615 | new_ref_elev = [] |
---|
| 616 | NODATA_value = -9999 |
---|
| 617 | for i in range(6): |
---|
| 618 | y = (6-i)*25.0 |
---|
| 619 | for j in range(5): |
---|
| 620 | x = j*25.0 |
---|
| 621 | z = x+2*y |
---|
| 622 | if j == 4: z = NODATA_value # column |
---|
| 623 | if i == 2 and j == 2: z = NODATA_value # random |
---|
| 624 | if i == 5 and j == 1: z = NODATA_value |
---|
| 625 | if i == 1: z = NODATA_value # row |
---|
| 626 | if i == 3 and j == 1: z = NODATA_value # two pts/row |
---|
| 627 | if i == 3 and j == 3: z = NODATA_value |
---|
| 628 | |
---|
| 629 | |
---|
| 630 | if z <> NODATA_value: |
---|
| 631 | new_ref_elev.append(z) |
---|
| 632 | new_ref_pts.append( [x,y] ) |
---|
| 633 | |
---|
| 634 | ref_points.append( [x,y] ) |
---|
| 635 | ref_elevation.append(z) |
---|
| 636 | |
---|
| 637 | fid.write('%f ' %z) |
---|
| 638 | fid.write('\n') |
---|
| 639 | |
---|
| 640 | fid.close() |
---|
| 641 | |
---|
| 642 | |
---|
| 643 | #Write prj file with metadata |
---|
| 644 | metafilename = root+'.prj' |
---|
| 645 | fid = open(metafilename, 'w') |
---|
| 646 | |
---|
| 647 | |
---|
| 648 | fid.write("""Projection UTM |
---|
| 649 | Zone 56 |
---|
| 650 | Datum WGS84 |
---|
| 651 | Zunits NO |
---|
| 652 | Units METERS |
---|
| 653 | Spheroid WGS84 |
---|
| 654 | Xshift 0.0000000000 |
---|
| 655 | Yshift 10000000.0000000000 |
---|
| 656 | Parameters |
---|
| 657 | """) |
---|
| 658 | fid.close() |
---|
| 659 | |
---|
| 660 | #Convert to NetCDF pts |
---|
| 661 | convert_dem_from_ascii2netcdf(root) |
---|
| 662 | dem2pts(root) |
---|
| 663 | |
---|
| 664 | #Check contents |
---|
| 665 | #Get NetCDF |
---|
| 666 | fid = NetCDFFile(root+'.pts', 'r') |
---|
| 667 | |
---|
| 668 | # Get the variables |
---|
| 669 | #print fid.variables.keys() |
---|
| 670 | points = fid.variables['points'] |
---|
| 671 | elevation = fid.variables['elevation'] |
---|
| 672 | |
---|
| 673 | #Check values |
---|
| 674 | #print 'points', points[:] |
---|
| 675 | assert len(points) == len(new_ref_pts), 'length of returned points not correct' |
---|
| 676 | assert allclose(points, new_ref_pts), 'points do not align' |
---|
| 677 | |
---|
| 678 | #print 'elevation', elevation[:] |
---|
| 679 | assert len(elevation) == len(new_ref_elev), 'length of returned elevation not correct' |
---|
| 680 | assert allclose(elevation, new_ref_elev), 'elevations do not align' |
---|
| 681 | |
---|
| 682 | #Cleanup |
---|
| 683 | fid.close() |
---|
| 684 | |
---|
| 685 | |
---|
| 686 | os.remove(root + '.pts') |
---|
| 687 | os.remove(root + '.dem') |
---|
| 688 | os.remove(root + '.asc') |
---|
| 689 | os.remove(root + '.prj') |
---|
| 690 | |
---|
| 691 | def NOT_test_dem2pts_bounding_box_Nullvalues(self): |
---|
| 692 | """Test conversion from dem in ascii format to native NetCDF xya format |
---|
| 693 | """ |
---|
| 694 | |
---|
| 695 | import time, os |
---|
| 696 | from Numeric import array, zeros, allclose, Float, concatenate |
---|
| 697 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 698 | |
---|
| 699 | #Write test asc file |
---|
| 700 | root = 'demtest' |
---|
| 701 | |
---|
| 702 | filename = root+'.asc' |
---|
| 703 | fid = open(filename, 'w') |
---|
| 704 | fid.write("""ncols 5 |
---|
| 705 | nrows 6 |
---|
| 706 | xllcorner 2000.5 |
---|
| 707 | yllcorner 3000.5 |
---|
| 708 | cellsize 25 |
---|
| 709 | NODATA_value -9999 |
---|
| 710 | """) |
---|
| 711 | #Create linear function |
---|
| 712 | |
---|
| 713 | ref_points = [] |
---|
| 714 | ref_elevation = [] |
---|
| 715 | new_ref_pts1 = [] |
---|
| 716 | new_ref_elev1 = [] |
---|
| 717 | NODATA_value = -9999 |
---|
| 718 | for i in range(6): |
---|
| 719 | y = (6-i)*25.0 |
---|
| 720 | for j in range(5): |
---|
| 721 | x = j*25.0 |
---|
| 722 | z = x+2*y |
---|
| 723 | if j == 4: z = NODATA_value # column |
---|
| 724 | if i == 2 and j == 2: z = NODATA_value # random |
---|
| 725 | if i == 5 and j == 1: z = NODATA_value |
---|
| 726 | if i == 1: z = NODATA_value # row |
---|
| 727 | if i == 3 and j == 1: z = NODATA_value # two pts/row |
---|
| 728 | if i == 3 and j == 3: z = NODATA_value |
---|
| 729 | |
---|
| 730 | if z <> NODATA_value: |
---|
| 731 | new_ref_elev1.append(z) |
---|
| 732 | new_ref_pts1.append( [x,y] ) |
---|
| 733 | |
---|
| 734 | ref_points.append( [x,y] ) |
---|
| 735 | ref_elevation.append(z) |
---|
| 736 | fid.write('%f ' %z) |
---|
| 737 | fid.write('\n') |
---|
| 738 | |
---|
| 739 | fid.close() |
---|
| 740 | |
---|
| 741 | #Write prj file with metadata |
---|
| 742 | metafilename = root+'.prj' |
---|
| 743 | fid = open(metafilename, 'w') |
---|
| 744 | |
---|
| 745 | |
---|
| 746 | fid.write("""Projection UTM |
---|
| 747 | Zone 56 |
---|
| 748 | Datum WGS84 |
---|
| 749 | Zunits NO |
---|
| 750 | Units METERS |
---|
| 751 | Spheroid WGS84 |
---|
| 752 | Xshift 0.0000000000 |
---|
| 753 | Yshift 10000000.0000000000 |
---|
| 754 | Parameters |
---|
| 755 | """) |
---|
| 756 | fid.close() |
---|
| 757 | |
---|
| 758 | #Convert to NetCDF pts |
---|
| 759 | convert_dem_from_ascii2netcdf(root) |
---|
| 760 | dem2pts(root, easting_min=2010.0, easting_max=2110.0, |
---|
| 761 | northing_min=3035.0, northing_max=3125.5) |
---|
| 762 | |
---|
| 763 | #Check contents |
---|
| 764 | #Get NetCDF |
---|
| 765 | fid = NetCDFFile(root+'.pts', 'r') |
---|
| 766 | |
---|
| 767 | # Get the variables |
---|
| 768 | #print fid.variables.keys() |
---|
| 769 | points = fid.variables['points'] |
---|
| 770 | elevation = fid.variables['elevation'] |
---|
| 771 | |
---|
| 772 | #Check values |
---|
| 773 | assert fid.xllcorner[0] == 2010.0 |
---|
| 774 | assert fid.yllcorner[0] == 3035.0 |
---|
| 775 | |
---|
| 776 | #create new reference points |
---|
| 777 | ref_points = [] |
---|
| 778 | ref_elevation = [] |
---|
| 779 | new_ref_pts2 = [] |
---|
| 780 | new_ref_elev2 = [] |
---|
| 781 | for i in range(4): |
---|
| 782 | y = (4-i)*25.0 + 25.0 |
---|
| 783 | y_new = y + 3000.5 - 3035.0 |
---|
| 784 | for j in range(4): |
---|
| 785 | x = j*25.0 + 25.0 |
---|
| 786 | x_new = x + 2000.5 - 2010.0 |
---|
| 787 | z = x+2*y |
---|
| 788 | |
---|
| 789 | if j == 3: z = NODATA_value # column |
---|
| 790 | if i == 1 and j == 1: z = NODATA_value # random |
---|
| 791 | if i == 4 and j == 0: z = NODATA_value |
---|
| 792 | if i == 0: z = NODATA_value # row |
---|
| 793 | if i == 2 and j == 0: z = NODATA_value # two pts/row |
---|
| 794 | if i == 2 and j == 2: z = NODATA_value |
---|
| 795 | |
---|
| 796 | if z <> NODATA_value: |
---|
| 797 | new_ref_elev2.append(z) |
---|
| 798 | new_ref_pts2.append( [x_new,y_new] ) |
---|
| 799 | |
---|
| 800 | |
---|
| 801 | ref_points.append( [x_new,y_new] ) |
---|
| 802 | ref_elevation.append(z) |
---|
| 803 | |
---|
| 804 | #print points[:] |
---|
| 805 | #print ref_points |
---|
| 806 | #assert allclose(points, ref_points) |
---|
| 807 | |
---|
| 808 | #print attributes[:] |
---|
| 809 | #print ref_elevation |
---|
| 810 | #assert allclose(elevation, ref_elevation) |
---|
| 811 | |
---|
| 812 | |
---|
| 813 | assert len(points) == len(new_ref_pts2), 'length of returned points not correct' |
---|
| 814 | assert allclose(points, new_ref_pts2), 'points do not align' |
---|
| 815 | |
---|
| 816 | #print 'elevation', elevation[:] |
---|
| 817 | assert len(elevation) == len(new_ref_elev2), 'length of returned elevation not correct' |
---|
| 818 | assert allclose(elevation, new_ref_elev2), 'elevations do not align' |
---|
| 819 | #Cleanup |
---|
| 820 | fid.close() |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | os.remove(root + '.pts') |
---|
| 824 | os.remove(root + '.dem') |
---|
| 825 | os.remove(root + '.asc') |
---|
| 826 | os.remove(root + '.prj') |
---|
| 827 | |
---|
[2344] | 828 | #******************** |
---|
| 829 | #*** END OF OBSOLETE FUNCTIONS |
---|
| 830 | #*************** |
---|