dtm to heightmap converter

dtm_to_heightmap.py

#!/usr/bin/env python

import re, os
from struct import *

# http://wms.lroc.asu.edu/lroc/view_rdr/NAC_DTM_APOLLO15
FileName = 'NAC_DTM_APOLLO15_E261N0036.IMG'
f = open(FileName, 'rb')

# reading meta info
line="BEGIN"
while line != "END":
    line = f.readline().replace('\r\n', '')
    if line.find('MAP_RESOLUTION') > 0:
        MapResolution = float(re.search(r"MAP_RESOLUTION\s*=\s*([\d\.]*)\s*", line).group(1))
    if line.find('MAP_SCALE') > 0:
        MapScale = float(re.search(r"MAP_SCALE\s*=\s*([\d\.]*)\s*", line).group(1))
    if line.find('MAXIMUM_LATITUDE') > 0:
        MaximumLatitude = float(re.search(r"MAXIMUM_LATITUDE\s*=\s*([\d\.]*)\s*", line).group(1))
    if line.find('MINIMUM_LATITUDE') > 0:
        MinimumLatitude = float(re.search(r"MINIMUM_LATITUDE\s*=\s*([\d\.]*)\s*", line).group(1))
    if line.find('EASTERNMOST_LONGITUDE') > 0:
        EasternmostLongitude = float(re.search(r"EASTERNMOST_LONGITUDE\s*=\s*([\d\.]*)\s*", line).group(1))
    if line.find('WESTERNMOST_LONGITUDE') > 0:
        WesternmostLongitude = float(re.search(r"WESTERNMOST_LONGITUDE\s*=\s*([\d\.]*)\s*", line).group(1))
    print line


print("--- Parsed Data")
print("Map Resolution {} <PIX/DEG>".format(MapResolution))
print("Map Sale {} <METERS>".format(MapScale))
print("Maxmum Latitude {} <DEG>".format(MaximumLatitude))
print("Minimum Latitude {} <DEG>".format(MinimumLatitude))
print("Easternmost Longitude {} <DEG>".format(EasternmostLongitude))
print("Westernmost Longitude {} <DEG>".format(WesternmostLongitude))

LinePixel= int(round((EasternmostLongitude - WesternmostLongitude)*MapResolution))
SamplePixel = int(round((MaximumLatitude - MinimumLatitude)*MapResolution))
print("... {} pixel / line".format(LinePixel))
print("... {} line / sample".format(SamplePixel))
print "-- Done reading headers"
    
while f.read(1) == ' ':
    pass
f.seek(-1,1)# unread
#f.readline()

# CORE_NULL                  = 16#FF7FFFFB#
# CORE_LOW_REPR_SATURATION   = 16#FF7FFFFC#
# CORE_LOW_INSTR_SATURATION  = 16#FF7FFFFD#
# CORE_HIGH_REPR_SATURATION  = 16#FF7FFFFF#
# CORE_HIGH_INSTR_SATURATION = 16#FF7FFFFE#
SkipData = (pack('BBBB', 0xFB, 0xFF, 0x7F, 0xFF), pack('BBBB', 0xFC, 0xFF, 0x7F, 0xFF), pack('BBBB', 0xFD, 0xFF, 0x7F, 0xFF), pack('BBBB', 0xFF, 0xFF, 0x7F, 0xFF), pack('BBBB', 0xFE, 0xFF, 0x7F, 0xFF))
altitude = 0

HeightMapSize = 257 # 257 pixel x 257 pixel (must be 2^n+1)
HeightMap = [[0 for i in range(HeightMapSize)] for j in range(HeightMapSize)]

# https://www.hq.nasa.gov/alsj/alsjcoords.html
CenterLatitude = 26.13224
CenterLongitude = 3.63400
# http://tools.wmflabs.org/geohack/geohack.php?pagename=Apollo_15&params=26_7_55.99_N_3_38_1.90_E_type:landmark_globe:moon&title=Apollo+15+landing
CenterLatitude = 26.132219 - 0.002
CenterLongitude = 3.633861
# http://www.hq.nasa.gov/alsj/a15/images15.html#11450
# http://commons.wikimedia.org/wiki/File:Apollo_15_Station_2_Rille,_Lunar_Rover,_Scott.jpg
# http://www.google.com/moon/#lat=26.103944&lon=3.568008&zoom=15&apollo=a15/26
CenterLatitude = 26.103944
CenterLongitude = 3.568008
# Spur Crater
# http://www.google.com/moon/#lat=25.962926&lon=3.675353&zoom=15&apollo=a15/18
# http://the-moon.wikispaces.com/Spur
CenterLatitude = 25.95
CenterLongitude = 3.624
# elbow
# http://the-moon.wikispaces.com/Elbow
# http://www.google.com/moon/#lat=26.008005&lon=3.604091&zoom=15&apollo=a15/6
# 1488 / 7622 pixel
CenterLatitude = 26.007005 + 0.025
CenterLongitude = 3.604091 - 0.00
# Dune Crater (https://www.hq.nasa.gov/alsj/a15/a15.trvsta4.html)
# https://www.hq.nasa.gov/alsj/a15/a15pan1463033s.jpg
# http://www.google.com/moon/#lat=26.006452&lon=3.666045&zoom=15&apollo=a15/20
# 2320 7588
CenterLatitude = 26.006452 + 0.0284
CenterLongitude = 3.666045 - 0.013

print("Applo 15 landing locaiton {} / {}".format(CenterLatitude, CenterLongitude))
StartLatitude = CenterLatitude + HeightMapSize/2/MapResolution
StartLongitude =  CenterLongitude - HeightMapSize/2/MapResolution
EndLatitude = CenterLatitude - HeightMapSize/2/MapResolution
EndLongitude =  CenterLongitude + HeightMapSize/2/MapResolution
print("HightMap Cropping Area {} - {}".format(StartLatitude, EndLatitude))
print("                       {} - {}".format(StartLongitude, EndLongitude))
if StartLatitude > MaximumLatitude or StartLatitude < MinimumLatitude or StartLongitude < WesternmostLongitude or EndLongitude > EasternmostLongitude:
    print("ERROR: Cropping  outside of the map")
    exit(1)

OffsetLatitude = int(round((MaximumLatitude - StartLatitude)*MapResolution))
OffsetLongitude = int(round((StartLongitude - WesternmostLongitude)*MapResolution))

print("HightMap Cropping Pixel Latitude  {} - {} - {}".format(OffsetLatitude, OffsetLatitude+HeightMapSize/2, OffsetLatitude+HeightMapSize))
print("                        Longitude {} - {} - {}".format(OffsetLongitude, OffsetLongitude+HeightMapSize/2, OffsetLongitude+HeightMapSize) )

MinimumAltitude =  20000
MaximumAltitude = -20000
for latitude in range(SamplePixel):
    for longitude in range(LinePixel):
        data = f.read(4)
        if len(data) != 4:
            print("could not read correct data ...", latitude, longitude)
        if data in SkipData:
            continue
        altitude = unpack('f', data)[0]
        if OffsetLatitude <= latitude and latitude < OffsetLatitude + HeightMapSize and \
           OffsetLongitude <= longitude and longitude < OffsetLongitude + HeightMapSize:
            HeightMap[latitude-OffsetLatitude][longitude-OffsetLongitude] = altitude
            # max/min withih HeightMap
            if altitude < MinimumAltitude:
                MinimumAltitude = altitude
            if altitude > MaximumAltitude:
                MaximumAltitude = altitude
print "-- Done reading data"
print len(f.read()), "data remaining ...."
print "Minimum Altitude MinimumAltitude", MinimumAltitude
print "Maximul Altitude MaximumAltitude", MaximumAltitude
HeightResolution = (MaximumAltitude - MinimumAltitude)/255
print "Altitude resolution ", HeightResolution

print "HeightMap Size ", HeightMapSize*MapScale, MaximumAltitude - MinimumAltitude
print "HeightMap Center Altitude ", (HeightMap[HeightMapSize/2][HeightMapSize/2]-MinimumAltitude)

for j in range(HeightMapSize):
    for i in range(HeightMapSize):
        if HeightMap[j][i] != 0:
            if int((HeightMap[j][i]-MinimumAltitude)/HeightResolution) < 0:
                print j, i, HeightMap[j][i], MinimumAltitude, HeightResolution, (HeightMap[j][i]-MinimumAltitude)/HeightResolution
            if int((HeightMap[j][i]-MinimumAltitude)/HeightResolution) > 255:
                print j, i, HeightMap[j][i], MinimumAltitude, HeightResolution, (HeightMap[j][i]-MinimumAltitude)/HeightResolution
        if HeightMap[j][i] != 0:
            HeightMap[j][i] = int((HeightMap[j][i]-MinimumAltitude)/HeightResolution)

OutputName = "{}_{}x{}+{}+{}.png".format(os.path.splitext(FileName)[0], HeightMapSize, HeightMapSize, OffsetLatitude, OffsetLongitude)
print  "-- Output", OutputName

# curl -LO https://raw.github.com/drj11/pypng/master/code/png.py
import png
f = open(OutputName, 'wb')
w = png.Writer(HeightMapSize, HeightMapSize, greyscale=True)
w.write(f, HeightMap)
f.close()
print "-- Done writing"