473,320 Members | 1,870 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,320 software developers and data experts.

quick script to read digital terrain elevation data?


I'm running python 2.3 on Windows XP. Anyone have a quick small script
to convert .DT1 and .DEM data to ASCII or some other format? I don't
need a viewer.

Thanks!

Feb 22 '06 #1
1 3961
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
elevationPattern = None
groundSystem = None
groundZone = None
projectParams = []
planeUnitOfMeasure = None
elevUnitOfMeasure = None
polygonSizes = None
groundCoords = {}
elevBounds = None
localRotation = None
accuracyCode = None

#((),(),())
spatialResolution = ()

profileDimension = None
unknownA = None
unknownB = None
unknownC = None
unknownD = None
unknownE = None
unknownF = None
unknownG = None
unknownH = None

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

self.cellCount = 0

for r in range(self.firstRow, self.lastRow, self.rowInt):
#scale the raw value
tempFloat = float(self.base[r]) * self.verticalScale
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(self, 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.elevationPattern = int(fh.read(6))
print "elevation Pattern: %i" % self.elevationPattern

self.groundSystem = int(fh.read(6))
print "planimetric reference system: %i" % self.groundSystem

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.projectParams.append(float(value.replace("D",
"E")))
print "proj: %15.7f %15.7f %15.7f" % (self.projectParams[-3],
self.projectParams[-2], self.projectParams[-1])

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

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

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

#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(24).replace("D", "E"))
lng = float(fh.read(24).replace("D", "E"))

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

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

self.localRotation = float(fh.read(24).replace("D", "E"))

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

self.spatialResolution = (
float(fh.read(12)),
float(fh.read(12)),
float(fh.read(12))
)

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

self.spatialResolution[1],

self.spatialResolution[2]
)

self.profileDimension = (
int(fh.read(6)),
int(fh.read(6))
)
print "map size is %i x %i" % (self.profileDimension[0],
self.profileDimension[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.columnCount):
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(profileSize)
#profileSize = (int(fh.read(6)), int(fh.read(6)))

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

localElevation = float(fh.read(24).replace("D", "E"))
elevExtremea = (
float(fh.read(24).replace("D", "E")),
float(fh.read(24).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.southMostSample) /
self.spatialResolution[1]

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

#skip ahead to the first row
self.base = []
for r in range(self.firstRow):
self.base.append(0)

for r in range(self.firstRow, self.lastRow):
value = fh.read(6)
print self.firstRow, r, self.lastRow, repr(value)
try:
self.base.append(int(value))
except:
raise "the horrors of war!"
#self.base.append(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(oh, 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.spatialResolution[0] == 3.0:
comp = 2.0
D.deltaY = D.spatialResolution[0] / 3600
else:
comp = 0.0
D.deltaY = D.spatialResolution[0]

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

#trunc' up to nearest spatialResolution
D.eastMostSample = int(D.eastMost / D.spatialResolution[0]) *
int(D.spatialResolution[0])
D.westMostSample = int((D.westMost + comp) / D.spatialResolution[0]) *
int(D.spatialResolution[0])
D.northMostSample = int(D.northMost / D.spatialResolution[1]) *
int(D.spatialResolution[1])
D.southMostSample = int((D.southMost + comp) / D.spatialResolution[1])
* int(D.spatialResolution[1])

if D.spatialResolution == 3.0:
#it's a 1x1 degree DEM
print "eastMostSample: %10f, westMostSample: %10f" %
(D.eastMostSample / 3600, D.westMostSample / 3600)
print "northMostSample: %10f, southMostSample: %10f" %
(D.northMostSample / 3600, D.southMostSample / 3600)

D.columnCount = (D.eastMostSample - D.westMostSample) /
int(D.spatialResolution[0]) + 1
D.rowCount = (D.northMostSample - D.southMostSample) /
int(D.spatialResolution[1]) + 1

if D.columnCount != D.profileDimension[1]:
print "CALCULATED column Count %i != HEADER %i" % (D.columnCount,
D.profileDimension[1])
print "will use SMALLER column Count"
D.columnCount = min(D.profileDimension[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.xyz"
oh = open(outfile, "w")

D.processProfiles(fh, oh)

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

Feb 23 '06 #2

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

0
by: DolphinDB | last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation. Take...
0
by: DolphinDB | last post by:
Tired of spending countless mintues downsampling your data? Look no further! In this article, you’ll learn how to efficiently downsample 6.48 billion high-frequency records to 61 million...
0
by: ryjfgjl | last post by:
ExcelToDatabase: batch import excel into database automatically...
0
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
1
isladogs
by: isladogs | last post by:
The next Access Europe meeting will be on Wednesday 6 Mar 2024 starting at 18:00 UK time (6PM UTC) and finishing at about 19:15 (7.15PM). In this month's session, we are pleased to welcome back...
0
by: Vimpel783 | last post by:
Hello! Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
0
by: Defcon1945 | last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
0
by: af34tf | last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you
0
by: Faith0G | last post by:
I am starting a new it consulting business and it's been a while since I setup a new website. Is wordpress still the best web based software for hosting a 5 page website? The webpages will be...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.