What’s the Point? Interpolation & Extrapolation with a Regular Grid DEM

1. What’s the Point?Interpolation & Extrapolation with a Regular Grid DEM David Kidner, Mark Dorey & Derek Smith University of Glamorgan School of Computing Pontypridd WALES, U.K. CF37 1DL e-mail: dbkidner@glam.ac.uk

2. What’s the Point? • Digital Terrain Modelling and “Grids” • What’s the Point ? • Interpolation Algorithms • Tests and Results • Extrapolation Algorithms • Data Compression • Tests and Results • Conclusions

3. A Digital Elevation Model (DEM) • Regular grid of elevations represents the heights at discrete samples of a continuous surface • vertices are sampled or interpolated independently • represented in a 2D matrix • No direct topological relationship between points • 2D Grid imposes an implicit representation of surface form • usually as a linear relationship between vertices • Simple and convenient

4. Which One’s the DEM? (a) Discrete Elevation Samples (b) Implicit (Linear) Continuous Surface

5. Interpolation • DEM resolution should be dependent upon the variability of each terrain surface • … but rarely is • The requirement of the DEM is to represent the terrain surface such that elevations can be retrieved or inferred for any given location • i.e. usually by interpolation • The method of interpolation is often ignored • Required for most, if not all applications

6. 68 54 D 50 72 What’s the Point? • Does Interpolation matter? • What’s the height at D ? • Where’s the 60m Contour(s) ?

7. 68 54 0.4 B D C 0.6 A 50 72 0.2 0.8 Interpolation for Visibility Analysis • What’s the profile through the “cell” ?

8. Interpolation for Visibility Analysis (a) Linear with Diagonal (b) Linear without Diagonal (c) Bilinear 20m Object (a) Completely Obscured (b) Completely Visible (c) Partially Visible

9. Interpolation Algorithms • Very small interpolation errors may lead to very large application errors • visibility analysis, hydrological modelling, contouring, etc. • Interpolation is flawed if we only consider the grid cell of the point to be estimated • Most GIS only consider the 4 vertices of the grid cell ! • bilinear interpolation

10. Interpolation Alternatives • For the most part, we can use polynomial interpolation of the form: hi = a00 + a10x + a01y + a20x2 + a11xy + a02y2 + a30x3 + a21x2y + a12xy2 + a03y3 + a31x3y + a22x2y2 + a13xy3 + a32x3y2 + a23x2y3 + a33x3y3 + … + amnxmyn • solved from the set of simultaneous equations which are set up, one for each point.

11. Interpolation Alternatives • Level Plane (1 coefficient) • Linear Plane (3) • Double Linear and Bilinear (4) • Biquadratic (8 or 9) • Bicubic (12 or 16) • Biquintic (36) • Jancaitis Biquadratic, Piecewise Cubics, etc.

12. 288 278 253 252 264 250 230 264 232 232 249 277 218 233 267 287 Linear 1 Linear 2 Double Linear Bilinear Biquadratic Bicubic Jancaitis Biquintic (9 term) (16 term) (36 term)

13. Results (1) Test Surface Functions (Franke, 1979; Akima, 1996)

14. Results (1) Test Surface Functions

15. Results (1) Test Surface Functions

16. Results (1) Test Surface Functions • Higher-order interpolation algorithms will always out-perform linear techniques

17. Results (2) Actual Terrain Based on Ordnance Survey data for S. Wales: • 1:50,000 Scale (50 m) DEMs and 1:10,000 Scale (10 m) DEMs • Higher-order interpolation algorithms will always out-perform linear techniques • By 3% to 10% (of the RMSE) • Less correlation as to which algorithm performs best

18. Extrapolation • Interpolation outside the spatial extent • Extrapolation can be considered to be at the heart of the best techniques for spatial data compression • i.e. what is the next symbol in the series • or “standing on the surface and given my field of view, what is the elevation if I take one step backwards?”

19. O.S. 50m DEM O.S. 10m DEM Why do we need DEM compression? • Seen as yesterday’s problem ? • expensive hardware; small capacity disks, etc. • File/Internet Transfer • Higher Resolutions 2m LiDAR DEM

20. DEM Extrapolation & Prediction

21. ORIGINAL DEM 10% Frequency 5% Elevation Range DEM OF ERROR CORRECTIONS 40% Frequency 20% Elevation Correction Range (Prediction Errors) DEM Transformation for Compression

22. TerrainExtrapolators

23. Linear Predictor: Z = a + g - f 343 338 333 332 327 337 343 347 352 358 369 375 381 393 405 413 420 420 421 356 352 349 345 339 350 351 361 359 360 365 373 378 386 397 406 413 417 418 378 368 364 359 354 360 364 369 367 364 369 372 378 385 393 401 408 413 413 393 382 375 368 363 369 372 375 377 373 371 375 379 384 390 397 402 409 409 403 394 386 378 373 375 378 381 384 385 380 379 382 385 390 394 399 402 403 413 407 396 388 381 381 384 386 389 391 390 384 390 390 390 393 395 397 398 423 418 409 399 392 390 390 391 392 392 391 388 390 390 391 391 392 394 395 432 426 420 409 403 399 398 398 396 395 393 390 390 390 391 390 389 390 391 a Z 438 435 428 422 414 410 405 402 399 398 395 391 389 387 385 386 385 385 384 450 448 443 433 425 419 410 407 402 401 398 393 388 383 380 378 379 379 377 455 454 453 447 441 432 426 420 409 405 401 395 388 382 377 370 372 373 370 460 457 457 454 450 444 438 430 419 408 404 398 390 384 378 370 366 368 364 460 459 460 459 458 451 448 442 431 417 408 402 393 387 381 372 366 361 359 458 459 461 463 462 459 456 452 442 427 414 405 396 388 382 374 367 361 357 f g 452 458 462 465 468 465 461 457 449 438 421 409 397 389 382 375 368 363 358 448 454 461 467 470 470 464 461 454 445 430 413 400 389 383 375 369 364 359 442 449 456 462 469 470 466 461 457 450 440 420 406 390 382 376 370 366 361 e.g. 462 = Z = 458 + 461 - 454 = 465 (Correction = -3) 343 338 -2 3 1 -1 5 -6 7 5 6 -2 1 4 1 -1 0 420 421 356 352 1 1 -1 5 -3 5 0 4 0 5 -1 1 3 1 0 417 418 378 368 3 2 0 0 1 2 -4 1 7 -1 2 2 2 1 2 413 413 393 382 1 1 0 4 0 0 -1 -5 3 5 1 2 1 3 0 409 409 403 394 3 0 2 2 0 1 0 -1 -4 5 -3 3 5 1 3 402 403 413 407 -2 2 0 2 3 1 2 2 0 -3 4 0 -1 3 1 397 398 423 418 -3 1 -1 2 1 1 3 1 1 0 2 0 0 1 2 394 395 432 426 1 -5 2 0 4 3 1 0 1 1 2 2 3 -2 0 390 391 438 435 -2 4 0 2 4 0 2 0 0 1 3 3 1 3 -2 385 384 450 448 -4 -4 -2 3 -3 3 6 3 1 1 2 1 2 5 -1 379 377 455 454 -1 -3 -2 -3 0 2 0 7 0 0 1 0 1 1 6 373 370 460 457 -1 -2 -3 1 -3 -2 0 3 5 0 1 0 0 1 2 368 364 460 459 -1 -3 0 -4 0 -2 -1 1 4 3 0 2 0 -1 1 361 359 458 459 -2 -1 -4 0 1 0 -2 -4 4 3 3 0 1 -1 0 361 357 452 458 -3 -3 0 -3 2 -1 -1 -2 -2 5 1 3 -1 1 -1 363 358 448 454 461 467 470 470 464 461 454 445 430 413 400 389 383 375 369 364 359 442 449 456 462 469 470 466 461 457 450 440 420 406 390 382 376 370 366 361

24. 343 338 333 332 327 337 343 347 352 358 369 375 381 393 405 413 420 420 421 356 352 349 345 339 350 351 361 359 360 365 373 378 386 397 406 413 417 418 12-Point Predictor: 378 368 364 359 354 360 364 369 367 364 369 372 378 385 393 401 408 413 413 393 382 375 368 363 369 372 375 377 373 371 375 379 384 390 397 402 409 409 403 394 386 378 373 375 378 381 384 385 380 379 382 385 390 394 399 402 403 413 407 396 388 381 381 384 386 389 391 390 384 390 390 390 393 395 397 398 423 418 409 399 392 390 390 391 392 392 391 388 390 390 391 391 392 394 395 b a Z 432 426 420 409 403 399 398 398 396 395 393 390 390 390 391 390 389 390 391 438 435 428 422 414 410 405 402 399 398 395 391 389 387 385 386 385 385 384 450 448 443 433 425 419 410 407 402 401 398 393 388 383 380 378 379 379 377 f g e h i 455 454 453 447 441 432 426 420 409 405 401 395 388 382 377 370 372 373 370 460 457 457 454 450 444 438 430 419 408 404 398 390 384 378 370 366 368 364 460 459 460 459 458 451 448 442 431 417 408 402 393 387 381 372 366 361 359 458 459 461 463 462 459 456 452 442 427 414 405 396 388 382 374 367 361 357 l m n o p 452 458 462 465 468 465 461 457 449 438 421 409 397 389 382 375 368 363 358 448 454 461 467 470 470 464 461 454 445 430 413 400 389 383 375 369 364 359 442 449 456 462 469 470 466 461 457 450 440 420 406 390 382 376 370 366 361 e.g. Z = w1*458+w2*454+...+w12*469 = 464 (Correction = -2) 343 338 -1 3 -3 2 0 -4 2 1 5 -1 0 3 2 0 1 420 421 356 352 -2 -1 -5 4 -5 5 -1 0 -1 3 -2 -1 1 0 0 417 418 378 368 2 1 -3 0 0 3 -1 0 4 -1 1 1 1 0 0 413 413 393 382 0 -1 -3 3 0 0 0 -3 1 3 1 0 -1 1 -1 409 409 403 394 2 0 0 1 -1 -1 -2 -2 -4 -1 -3 -1 1 0 2 402 403 413 407 -2 1 -1 1 2 0 1 1 0 -4 4 -1 -1 1 0 397 398 423 418 0 1 -1 0 -2 -1 1 0 0 -2 0 -3 -1 -1 0 394 395 432 426 2 -4 -1 -2 0 2 0 1 0 -1 0 -1 1 -2 -1 390 391 438 435 1 4 2 4 3 1 -1 0 -1 -1 1 0 -1 1 -1 385 384 450 448 -2 -3 -1 1 -4 0 0 1 1 1 1 1 2 2 0 379 377 455 454 2 -1 0 -2 1 2 -2 3 0 0 0 1 1 -2 2 373 370 460 457 1 0 0 1 0 -1 -1 -1 2 -1 -1 -1 0 -1 1 368 364 460 459 2 0 2 -4 0 0 0 -1 0 1 -2 1 1 -1 1 361 359 458 459 0 1 -2 -1 0 1 -1 -3 0 1 1 0 1 -1 0 361 357 452 458 -2 -2 1 -2 0 0 1 1 -1 2 -1 0 0 1 0 363 358 448 454 461 467 470 470 464 461 454 445 430 413 400 389 383 375 369 364 359 442 449 456 462 469 470 466 461 457 450 440 420 406 390 382 376 370 366 361

25. 60 50 F r 40 e q u 30 e n 20 c y 10 12-Point 4-Point 0 -17 -15 -14 -13 3-Point -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 Horizontal 0 1 2 3 4 5 6 Height Correction (m) 7 8 9 10 11 12 Frequency Distribution of 15x15 Corrections

26. O.S. South Wales 1201x801 DEM

27. Prediction (extrapolation) Corrections

29. Data Compression Results • GZIPped DEM requires a storage capacity of 261% of the best extrapolator and Arithmetic Coding method

30. Summary • Mathematical modelling has now largely been forgotten by today’s GIS developers • Many GIS techniques are of limited value and may propagate through to application error (e.g. visibility analysis) • For DEM Interpolation • don’t use linear algorithms • Mathematical modelling offers significant savings for spatial data compression