<commit_before>#!/usr/bin/python

from math import sqrt
import sys
import argparse
from collections import deque
from struct import pack, unpack

from osgeo import gdal

gdal.UseExceptions()
gdal.TermProgress = gdal.TermProgress_nocb

#
# NormalVector
#
# Calculate the normal vector of a triangle. (Unit vector perpendicular to
# triangle surface, pointing away from the "outer" face of the surface.)
#
# Parameters:
#  triangle vertices (nested x y z tuples)
#
# Returns:
#  normal vector (x y z tuple)
#
def NormalVector(t):
	
	(ax, ay, az) = t[0]
	(bx, by, bz) = t[1]
	(cx, cy, cz) = t[2]
	
	# first edge
	e1x = ax - bx
	e1y = ay - by
	e1z = az - bz
	
	# second edge
	e2x = bx - cx
	e2y = by - cy
	e2z = bz - cz
	
	# cross product
	cpx = e1y * e2z - e1z * e2y
	cpy = e1z * e2x - e1x * e2z
	cpz = e1x * e2y - e1y * e2x
	
	# return cross product vector normalized to unit length
	mag = sqrt((cpx * cpx) + (cpy * cpy) + (cpz * cpz))
	return (cpx/mag, cpy/mag, cpz/mag)

# stlwriter is a simple class for writing binary STL meshes.
# Class instances are constructed with a predicted face count.
# The output file header is overwritten upon completion with
# the actual face count.
class stlwriter():
	
	# path: output binary stl file path
	# facet_count: predicted number of facets
	def __init__(self, path, facet_count=0):
		
		self.f = open(path, 'wb')
		
		# track number of facets actually written
		self.written = 0
		
		# write binary stl header with predicted facet count
		self.f.write('\0' * 80)
		self.f.write(pack('<I', facet_count))
	
	# t: ((ax, ay, az), (bx, by, bz), (cx, cy, cz))
	def add_facet(self, t):
		# strictly speaking, we don't need to compute NormalVector,
		# as other tools could be used to update the output mesh.
		self.f.write(pack('<3f', *NormalVector(t)))
		for vertex in t:
			self.f.write(pack('<3f', *vertex))
		self.f.write('\0\0')
		self.written += 1
	
	def done(self):
		# update final facet count in header before closing file
		self.f.seek(80)
		self.f.write(pack('<I', self.written))
		self.f.close()
	
	def __enter__(self):
		return self
	
	def __exit__(self, exc_type, exc_value, traceback):
		self.done()

def fail(msg):
	print >> sys.stderr, msg
	exit(1)

def log(msg):
	if args.verbose:
		print >> sys.stderr, msg

ap = argparse.ArgumentParser(description='Convert a GDAL raster (like a GeoTIFF heightmap) to an STL terrain surface.')
ap.add_argument('-x', action='store', default=0.0, type=float, help='Fit output x to extent (mm)')
ap.add_argument('-y', action='store', default=0.0, type=float, help='Fit output y to extent (mm)')
optg_z = ap.add_mutually_exclusive_group()
optg_z.add_argument('-z', action='store', default=None, type=float, help='Z scale expressed as a vertical scale factor (1)')
optg_z.add_argument('-s', action='store', default=None, type=float, help='Z scale expressed as a ratio of vertical units per horizontal unit (1)')
ap.add_argument('-b', '--base', action='store', default=0.0, type=float, help='Base height (0)')
ap.add_argument('-c', '--clip', action='store_true', default=False, help='Clip z to minimum elevation')
ap.add_argument('-v', '--verbose', action='store_true', default=False, help='Print log messages')
ap.add_argument('--band', action='store', default=1, type=int, help='Raster data band (1)')
ap.add_argument('-m', '--minimum', action='store', default=None, type=float, help='Omit vertices below minimum elevation')
ap.add_argument('-M', '--maximum', action='store', default=None, type=float, help='Omit vertices above maximum elevation')
optg_region = ap.add_mutually_exclusive_group()
optg_region.add_argument('-w', '--window', action='store', default=None, type=float, nargs=4, help='Opposing corner coordinates in geographic CRS')
optg_region.add_argument('-p', '--pixels', action='store', default=None, type=float, nargs=4, help='Opposing corner coordinates in pixel coordinates')
ap.add_argument('RASTER', help='Input heightmap image')
ap.add_argument('STL',  help='Output STL path')
args = ap.parse_args()

try:
	img = gdal.Open(args.RASTER)
except RuntimeError, e:
	fail(str(e).strip())

# input raster dimensions
w = img.RasterXSize
h = img.RasterYSize
log("raster dimensions = (%s, %s)" % (str(w), str(h)))

# get default transformation from image coordinates to world coordinates
t = img.GetGeoTransform()

git = gdal.InvGeoTransform(t)[1]
log("geo->pixel transform: %s" % str(git))

# if a geographic window is specified, convert it to a pixel window in input raster coordinates
if args.window != None:
	
	# apply inverse geo transform to window points
	px0, py0 = gdal.ApplyGeoTransform(git, args.window[0], args.window[1])
	px1, py1 = gdal.ApplyGeoTransform(git, args.window[2], args.window[3])
	
	# set arg.pixels to obtained pixel points
	args.pixels = [px0, py0, px1, py1]

# if no extent window is specified, use full input raster.
if args.pixels == None:
	xmin = 0
	ymin = 0
	ww = w
	wh = h

# if a pixel extent window is specified (either explicitly or
# derived from a geographic window), clip it to the input raster extent.
else:
	
	xmin = int(round(min(args.pixels[0], args.pixels[2])))
	ymin = int(round(min(args.pixels[1], args.pixels[3])))
	
	xmax = int(round(max(args.pixels[0], args.pixels[2])))
	ymax = int(round(max(args.pixels[1], args.pixels[3])))
	
	if xmin >= w:
		fail("Region of interest lies entirely outside raster (xmin)")
	<commit_msg>Clarify input region window conditional statements<commit_after>#!/usr/bin/python

from math import sqrt
import sys
import argparse
from collections import deque
from struct import pack, unpack

from osgeo import gdal

gdal.UseExceptions()
gdal.TermProgress = gdal.TermProgress_nocb

#
# NormalVector
#
# Calculate the normal vector of a triangle. (Unit vector perpendicular to
# triangle surface, pointing away from the "outer" face of the surface.)
#
# Parameters:
#  triangle vertices (nested x y z tuples)
#
# Returns:
#  normal vector (x y z tuple)
#
def NormalVector(t):
	
	(ax, ay, az) = t[0]
	(bx, by, bz) = t[1]
	(cx, cy, cz) = t[2]
	
	# first edge
	e1x = ax - bx
	e1y = ay - by
	e1z = az - bz
	
	# second edge
	e2x = bx - cx
	e2y = by - cy
	e2z = bz - cz
	
	# cross product
	cpx = e1y * e2z - e1z * e2y
	cpy = e1z * e2x - e1x * e2z
	cpz = e1x * e2y - e1y * e2x
	
	# return cross product vector normalized to unit length
	mag = sqrt((cpx * cpx) + (cpy * cpy) + (cpz * cpz))
	return (cpx/mag, cpy/mag, cpz/mag)

# stlwriter is a simple class for writing binary STL meshes.
# Class instances are constructed with a predicted face count.
# The output file header is overwritten upon completion with
# the actual face count.
class stlwriter():
	
	# path: output binary stl file path
	# facet_count: predicted number of facets
	def __init__(self, path, facet_count=0):
		
		self.f = open(path, 'wb')
		
		# track number of facets actually written
		self.written = 0
		
		# write binary stl header with predicted facet count
		self.f.write('\0' * 80)
		self.f.write(pack('<I', facet_count))
	
	# t: ((ax, ay, az), (bx, by, bz), (cx, cy, cz))
	def add_facet(self, t):
		# strictly speaking, we don't need to compute NormalVector,
		# as other tools could be used to update the output mesh.
		self.f.write(pack('<3f', *NormalVector(t)))
		for vertex in t:
			self.f.write(pack('<3f', *vertex))
		self.f.write('\0\0')
		self.written += 1
	
	def done(self):
		# update final facet count in header before closing file
		self.f.seek(80)
		self.f.write(pack('<I', self.written))
		self.f.close()
	
	def __enter__(self):
		return self
	
	def __exit__(self, exc_type, exc_value, traceback):
		self.done()

def fail(msg):
	print >> sys.stderr, msg
	exit(1)

def log(msg):
	if args.verbose:
		print >> sys.stderr, msg

ap = argparse.ArgumentParser(description='Convert a GDAL raster (like a GeoTIFF heightmap) to an STL terrain surface.')
ap.add_argument('-x', action='store', default=0.0, type=float, help='Fit output x to extent (mm)')
ap.add_argument('-y', action='store', default=0.0, type=float, help='Fit output y to extent (mm)')
optg_z = ap.add_mutually_exclusive_group()
optg_z.add_argument('-z', action='store', default=None, type=float, help='Z scale expressed as a vertical scale factor (1)')
optg_z.add_argument('-s', action='store', default=None, type=float, help='Z scale expressed as a ratio of vertical units per horizontal unit (1)')
ap.add_argument('-b', '--base', action='store', default=0.0, type=float, help='Base height (0)')
ap.add_argument('-c', '--clip', action='store_true', default=False, help='Clip z to minimum elevation')
ap.add_argument('-v', '--verbose', action='store_true', default=False, help='Print log messages')
ap.add_argument('--band', action='store', default=1, type=int, help='Raster data band (1)')
ap.add_argument('-m', '--minimum', action='store', default=None, type=float, help='Omit vertices below minimum elevation')
ap.add_argument('-M', '--maximum', action='store', default=None, type=float, help='Omit vertices above maximum elevation')
optg_region = ap.add_mutually_exclusive_group()
optg_region.add_argument('-w', '--window', action='store', default=None, type=float, nargs=4, help='Opposing corner coordinates in geographic CRS')
optg_region.add_argument('-p', '--pixels', action='store', default=None, type=float, nargs=4, help='Opposing corner coordinates in pixel coordinates')
ap.add_argument('RASTER', help='Input heightmap image')
ap.add_argument('STL',  help='Output STL path')
args = ap.parse_args()

try:
	img = gdal.Open(args.RASTER)
except RuntimeError, e:
	fail(str(e).strip())

# input raster dimensions
w = img.RasterXSize
h = img.RasterYSize
log("raster dimensions = (%s, %s)" % (str(w), str(h)))

# get default transformation from image coordinates to world coordinates
t = img.GetGeoTransform()

git = gdal.InvGeoTransform(t)[1]
log("geo->pixel transform: %s" % str(git))

if args.window != None:
	# if a geographic window is specified, convert it to a pixel window in input raster coordinates
	
	# apply inverse geo transform to window points
	px0, py0 = gdal.ApplyGeoTransform(git, args.window[0], args.window[1])
	px1, py1 = gdal.ApplyGeoTransform(git, args.window[2], args.window[3])
	
	# set arg.pixels to obtained pixel points
	args.pixels = [px0, py0, px1, py1]

if args.pixels == None:
	# if no pixel extent window is specified, use whole input raster.
	xmin = 0
	ymin = 0
	ww = w
	wh = h
else:
	# if a pixel extent window is specified (either directly with --pixels or
	# derived from a geographic --window), clip it to the input raster extent.
	
	xmin = int(round(min(args.pixels[0], args.pixels[2])))
	ymin = int(round(min(args.pixels[1], args.pixels[3])))
	
	xmax = int(round(max(args.pixels[0], args.pixels[2])))
	ymax = int(round(max(args.pixels[1], args.pixels[3])))
	
	if xmin >= w:
		fail("Region of interest lies entirely outside raster (xmin)")
	
{'fraction_non_alphanumeric': 0.14318792257108456, 'fraction_numerical': 0.012966564786514773, 'mean_word_length': 4.526100307062436, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '840562', 'n_tokens_mistral': 3930, 'n_tokens_neox': 3735, 'n_words': 1194}