quick script to read digital terrain elevation data?

**Kane** · Feb 23 '06, 12:15 AM

Re: quick script to read digital terrain elevation data?

I have some in-house python for dem2xyz but it isn't pretty. I would
suggest starting with ftp.blm.gov/pub/gis/dem2xyz6.zip ; it is a C
program with source for just such conversions (my initial dem2xyz.py is
a direct python rewrite of the C code). I think it's actually easier
to read in C (when reading C I expect C conventions, when reading
python I expect Python, this is C written in python and it gives off an
ugly vibe)

Sorry so long, the task is more complicated than a quick small script.

Anyway, for what it is worth:
------------------------------------------------
#Python rewrite of 'DEM2XYZ.C' which isn't a friendly compile
# slower than C, sure. But it should work on Win and UNIX without any
hassles.
#if this reads like a C program you know why. Sorry ;)

####
## from 'DEM2XYZ6.ZIP pulled from ftp.blm.gov
#
# Origional Header comments:
# convert dem to x, y, z format, one dem profile at a time,
# converts 3 arc second to lat/long, sol katz, mar. 94,
# added sampling and cutting to size, sol katz, apr 94
# changed the calculation of start position in sampling code,
# sol katz, jan 95
# ver 4, got rid of last 160 bytes on header line. sol katz, mar
97

class DEM:
mapLabel = None
DEMlevel = None
elevationPatter n = None
groundSystem = None
groundZone = None
projectParams = []
planeUnitOfMeas ure = None
elevUnitOfMeasu re = None
polygonSizes = None
groundCoords = {}
elevBounds = None
localRotation = None
accuracyCode = None

#((),(),())
spatialResoluti on = ()

profileDimensio n = None
unknownA = None
unknownB = None
unknownC = None
unknownD = None
unknownE = None
unknownF = None
unknownG = None
unknownH = None

colStr = 0
def writeNormal(sel f, oh, XCoord, YCoord):
if self.spatialRes olution[1] == 3.0:
#3 degrees of arc
XCoord = XCoord / 3600.0
YCoord = YCoord / 3600.0

self.cellCount = 0

for r in range(self.firs tRow, self.lastRow, self.rowInt):
#scale the raw value
tempFloat = float(self.base[r]) * self.verticalSc ale
oh.write("%f %f %i\n" % (XCoord, YCoord, int(tempFloat)) )
self.cellCount += 1
#move up the delta y distance
YCoord += self.deltaY
return None

def writeSubset(sel f, oh, XCoord, YCoord):
return None

def getheader(self, fh):
print ""
self.mapLabel = fh.read(144)
print "Map label: %s" % self.mapLabel

self.DEMlevel = int(fh.read(6))
print "DEMlevel: %i" % self.DEMlevel

self.elevationP attern = int(fh.read(6))
print "elevation Pattern: %i" % self.elevationP attern

self.groundSyst em = int(fh.read(6))
print "planimetri c reference system: %i" % self.groundSyst em

self.groundZone = int(fh.read(6))
print "ground Zone: %i" % self.groundZone

#pull off 15, 24 byte sections and store in projectParams list
for k in range(5):
for l in range(3):
value = fh.read(24)
#we need to handle the exponent.. a simple float()
ain't doing the trick.
self.projectPar ams.append(floa t(value.replace ("D",
"E")))
print "proj: %15.7f %15.7f %15.7f" % (self.projectPa rams[-3],
self.projectPar ams[-2], self.projectPar ams[-1])

self.planeUnitO fMeasure = int(fh.read(6))
print "plane Unit Of Measure: %i" % self.planeUnitO fMeasure

self.elevUnitOf Measure = int(fh.read(6))
print "elevation measurement units code: %i" %
self.elevUnitOf Measure

self.polygonSiz es = int(fh.read(6))
print "polygon Sides: %i" % self.polygonSiz es

#nested list, ((lat,lng),(lat ,lng),(lat,lng) ,(lat,lng))
for k in ["NW", "NE", "SW", "SE"]:
#NW, NE, SW, SE
lat = float(fh.read(2 4).replace("D", "E"))
lng = float(fh.read(2 4).replace("D", "E"))

self.groundCoor ds[k] = (lat, lng)
print "%s (%15.5f, %15.5f)" % (k, lat, lng)

self.elevBounds = (
float(fh.read(2 4).replace("D", "E")),
float(fh.read(2 4).replace("D", "E"))
)

self.localRotat ion = float(fh.read(2 4).replace("D", "E"))

self.accuracyCo de = fh.read(6)
print "Min: %15.5f, Max: %15.5f, Accuracy Code: %s" %
(self.elevBound s[0], self.elevBounds[1], self.accuracyCo de)

self.spatialRes olution = (
float(fh.read(1 2)),
float(fh.read(1 2)),
float(fh.read(1 2))
)

print "Spatial Resolution: %15.5f, %15.5f, %15.5f" %
(self.spatialRe solution[0],

self.spatialRes olution[1],

self.spatialRes olution[2]
)

self.profileDim ension = (
int(fh.read(6)) ,
int(fh.read(6))
)
print "map size is %i x %i" % (self.profileDi mension[0],
self.profileDim ension[1])
#dump the next 160 bytes
self.unknownA = fh.read(6)
self.unknownB = fh.read(16)
self.unknownC = fh.read(2)
self.unknownD = fh.read(2)
self.unknownE = fh.read(2)
self.unknownF = fh.read(4)
self.unknownG = fh.read(4)
self.unknownH = fh.read(124)

return None

def processProfiles (self, fh, oh):
lastProfile = 0
for c in range(self.colu mnCount):
print "Working on column %s" % c
current, next = fh.read(6), fh.read(6)
try:
profileID = {"current": int(current), "next":
int(next)}
except:
raise "Failed: %s, %s" % (repr(current), repr(next))

print repr(profileID)
#profileID = (fh.read(6), int(fh.read(6)) )
profileSize = {"alpha": int(fh.read(6)) , "beta":
fh.read(6)}
print repr(profileSiz e)
#profileSize = (int(fh.read(6) ), int(fh.read(6)) )

planCoords = (
float(fh.read(2 4).replace("D", "E")),
float(fh.read(2 4).replace("D", "E")),
)
print "planCoords : %15.5f, %15.5f" % planCoords

localElevation = float(fh.read(2 4).replace("D", "E"))
elevExtremea = (
float(fh.read(2 4).replace("D", "E")),
float(fh.read(2 4).replace("D", "E"))
)

#'kludge to force the end of processing'???
#if (profileID["current"] - 1) != lastProfile:
# print "%s - 1 != %s" % (profileID["current"],
lastProfile)
# print "%i lines were written to the file" %
self.cellCount
# return None

lastProfile = profileID["next"]
print "Column %i has %i rows" % (profileID["next"],
profileSize["alpha"])
self.firstRow = int(planCoords[1] - self.southMostS ample) /
self.spatialRes olution[1]

self.lastRow = self.firstRow + profileSize["alpha"]

#skip ahead to the first row
self.base = []
for r in range(self.firs tRow):
self.base.appen d(0)

for r in range(self.firs tRow, self.lastRow):
value = fh.read(6)
print self.firstRow, r, self.lastRow, repr(value)
try:
self.base.appen d(int(value))
except:
raise "the horrors of war!"
#self.base.appe nd(int(fh.read( 3)))

#if cutting out a section, adjust the rows
if outType == 2:
#subset
self.firstRow = rowStr
planCoords[1] = planCoords[1] + (self.firstRow - 1) *
self.deltaY
lastRow = max(lastRow, rowEnd)

mod = c % colInt
if mod == 0 and c >= self.colStr:
if outType == 2:
writeSubset(oh, planCoords[0], planCoords[1])
else:
D.writeNormal(o h, planCoords[0], planCoords[1])

#trailer bytes?
#fh.read(424)

return None

YMAX = 2048
XMAX = 2048
SW = 0
NW = 1
NE = 2
SE = 3

infile = "/home/jkane/pypov-0.0.2/rawdata/1843.CDO"

print "DEM to x,y,z ascii file, Python variation based on C version 6"
print " Python version by Jason Kane, BroadLink Communications 2004"
print " C version by Sol Katz, BLM April 1997"

fh = open(infile)
D = DEM()
D.getheader(fh)

D.verticalScale = 1.0

#easy enough...
if D.spatialResolu tion[0] == 3.0:
comp = 2.0
D.deltaY = D.spatialResolu tion[0] / 3600
else:
comp = 0.0
D.deltaY = D.spatialResolu tion[0]

D.eastMost = max(D.groundCoo rds["NE"][0], D.groundCoords["SE"][0])
D.westMost = min(D.groundCoo rds["NW"][0], D.groundCoords["SW"][0])
D.northMost = max(D.groundCoo rds["NE"][1], D.groundCoords["NW"][1])
D.southMost = min(D.groundCoo rds["SW"][1], D.groundCoords["SE"][1])

#trunc' up to nearest spatialResoluti on
D.eastMostSampl e = int(D.eastMost / D.spatialResolu tion[0]) *
int(D.spatialRe solution[0])
D.westMostSampl e = int((D.westMost + comp) / D.spatialResolu tion[0]) *
int(D.spatialRe solution[0])
D.northMostSamp le = int(D.northMost / D.spatialResolu tion[1]) *
int(D.spatialRe solution[1])
D.southMostSamp le = int((D.southMos t + comp) / D.spatialResolu tion[1])
* int(D.spatialRe solution[1])

if D.spatialResolu tion == 3.0:
#it's a 1x1 degree DEM
print "eastMostSample : %10f, westMostSample: %10f" %
(D.eastMostSamp le / 3600, D.westMostSampl e / 3600)
print "northMostSampl e: %10f, southMostSample : %10f" %
(D.northMostSam ple / 3600, D.southMostSamp le / 3600)

D.columnCount = (D.eastMostSamp le - D.westMostSampl e) /
int(D.spatialRe solution[0]) + 1
D.rowCount = (D.northMostSam ple - D.southMostSamp le) /
int(D.spatialRe solution[1]) + 1

if D.columnCount != D.profileDimens ion[1]:
print "CALCULATED column Count %i != HEADER %i" % (D.columnCount,
D.profileDimens ion[1])
print "will use SMALLER column Count"
D.columnCount = min(D.profileDi mension[1], D.columnCount)

print "number of rows %i, number of columns %i" % (D.rowCount,
D.columnCount)
colInt = 1
D.rowInt = 1
outType = input("Enter 0 for all, 1 for samples, 2 for subset : ")

if outType == 1:
#samples
colInt = input("Enter Column sample interval : ")
D.rowInt = input("Enter row sample interval : ")
D.deltaY = D.deltaY * D.rowInt

elif outType == 2:
#subset
D.colStr = input("Enter start column (from left) : ")
colEnd = input("Enter end column (from left) : ")
rowStr = input("Enter start row (from bottom) : ")
rowEnd = input("Enter end row (from bottom) : ")

D.verticalScale = input("Enter a vertical scaling factor: ")

if (outType == 2 and D.columnCount > colEnd):
D.columnCount = colEnd

outfile = "outfile.xy z"
oh = open(outfile, "w")

D.processProfil es(fh, oh)

print "%i lines were written to the file." % D.cellCount

quick script to read digital terrain elevation data?

quick script to read digital terrain elevation data?

Comment