Displaying geographical elevation data

nomainwin

open "W020N90.DEM" for input as #fileIn
    dataFile$ =input$( #fileIn, lof( #fileIn))
    timer 20000, [o]
    wait
    [o]
    timer 0
close #fileIn

UpperLeftX   =  20
UpperLeftY   =  20
WindowWidth  =1250
WindowHeight =1050

open "Display data" for graphics_nf_nsb as #w

#w "trapclose [quit]"
#w "down ; fill black"

index =0
c     =0

'   data is HiLo. Elevations in metres 1 --> 9000 or so.
'   $d8 $f1 everywhere at sea level ( or no data available) ( doc'n says -9999????)
'   16-bit signed integer data in a simple binary raster, ie two bytes per coordinate.
'       $00 $03 represents +3m asl; $00 $23 represents +35m; $01 $23 represents +288m.
'   NROWS         6000
'   NCOLS         4800

'       The data are stored in row major order (all the data for row 1, followed by all the data for row 2, etc.).

ww =1700

x1 =1100: x2 =x1 +ww  '   select square area of UK
y1 =3450: y2 =y1 +ww

for downCol =y1 to y2 step ( y2 -y1) /1000
    for alongRow =x1 to x2 step ( x2 -x1) /800
        index     = 1 +int( alongRow) *2 +int( downCol) *9600
        highByte  = asc( mid$( dataFile$, index,    1))
        if highByte <128 then
            lowByte   = asc( mid$( dataFile$, index +1, 1))
            wordValue = 256 *highByte +lowByte
            'print alongRow, downCol, wordValue

            c =int( wordValue /4)
            col$ =str$( c); " 100 "; str$( 255 -c)
        else
            col$ ="170 220 255"
        end if
        #w "color "; col$
        #w "set "; 20 +800 *( alongRow -x1) /( x2 -x1); " "; 20 +1000 *( downCol -y1) /( y2 -y1)
        scan
    next alongRow
    #w "flush"
next downCol

#w "flush"

wait

[quit]
close #w
end