#!/usr/bin/python

#ddfex2.py - DDF Example #2 1.0 - 9 Feb 01
#Basic cross-platform script to extract SDTS DDF DEM data.
#See http://www.3dartist.com/WP/sdts/sdtsnotes.htm#datar
#  for an explanation.
#Update at http://www.3dartist.com/WP/python/pycode.htm#ddfex
#Created by Bill Allen <dempy@3dartist.com>

#HISTORY:
#2001-02-09 [1.0] First posted.

#PURPOSE:
#This script demos a basic jury-rigged approach to reading DEM
#  data from many 30 or 10m USGS SDTS DEM DDFs. It does provide
#  useful output in the way of a .raw 16-bit grayscale image
#  file that Photoshop can read.

#PLATFORM: This code should run on any machine that has Python
#  2.0 installed, but has been tested only under Windows ME.
#  Visit http://www.python.org to learn more about Python, a
#  cross-platform freeware language.

#LICENSE: This code is free software provided as is and without
#  warranty for fitness for any purpose. Use it at your own risk.
#  This code and derivations may be used anywhere, including for
#  commercial application. If used for public purposes, an
#  acknowledgment and a link back would be appreciated.

### USER-ASSIGNED VARIABLES

fn = 'c:/python/demtest/1168CEL0.DDF'   #put your DDF file's name here
ticker = 1                              #1 to report status, 0 for not

### BEGIN PROGRAM ###

print '\nReading elevations:',fn
f = open(fn,'rb')        #contains binary data, so read binary
r = f.read(7)
if r[6] != 'L' or fn[-8:].upper() != 'CEL0.DDF':
    print '\tNot a DDF file! ID =',r[6],'(\\',ord(r[6]),')\n'
else:                               #this appears to be a DDR
    r = r + f.read(int(r[0:5])-7)   #  so check the rest of it
    if ord(r[-1]) != 30:
        print '\tFailure reading DDR terminator:',ord(r[-1])
    elif r.find('(B(16))') < 0:
        print '\tUnknown data type:',r[-8:-1]
    else:                #we've got a good DDF, now read the records
        rc = 0
        if ticker:
            t = ''
            ts = '>'
        rFlag = 0
        g = open(fn[:-4]+'.raw','wb')  #file to receive elevations
        while 1:
            r = f.read(24)             #bring in DR leader 
            if not r:                  #end of file
                break
            rc = rc + 1
            if r[6] != 'D':
                if r[6] == 'R':
                    print '\tCan\'t handle leaderless records.\n'
                else:
                    print '\tNot a data record! ID =',r[6],'(',ord(r[6]),')\n'
                break    #bail out
            if rc == 1:                #do only once
                aPtr = int(r[12:17])   #  find data area starting point
                fLen = int(r[21]) + 1  #  position of # giving jump past aPtr 
            r = r + f.read(aPtr-24)    #read in the DR directory
            dPtr = int(r[aPtr-fLen:aPtr-1])  #calc jump past aPtr to the data
            r = r + f.read(dPtr)       #read to the end of attributes
            if rc == 1:
                roWidth = (int(r[0:5])-(aPtr + dPtr))-1   #width of data area
                colCnt = roWidth / 2
            d = f.read(roWidth)
            g.write(d)
            d = f.read(1)              #get the last byte of the record
            if ord(d) != 30:
                print '\tFailure reading DR terminator:',ord(d)
                break
            elif ticker:
                t = t + ts             #simple busy signal
                if len(t) == 79:
                    print t
                    if ts == '>':
                        ts = '<'
                    else:
                        ts = '>'
                    t = ''
        g.flush()
        g.close()
        print '\tCols =',colCnt,'Rows =',rc,'in 16-bit 1-channel \"Mac\" format'
f.close()