#! /usr/bin/python

# point_densometer.py   V 0.1
# copyright by Florian Lang under the GPL (www.gnu.org/licence)
#
# It is a "polished" piece of throw away code. For more info you might want
# to look on "http://home.arcor.de/florianlang/".
#
# You need the Python Image Libraries (PIL) - google for it.
#
# Usage: point_densometer.py FILE  ...where FILE is a textfile with coordinates 
# 
# Further parameters have to be set by adjusting the next three variables.
#
# The coordinates in the text file need to have the following format: 
# 1000    345235.56  234513.64  345.32
# number  x          y          alt
# 
# That will work for sure, however you can use any format - just adopt some lines
# in this file - search for FNORD to find the place and a description.

gridspacing = 5                 # the side length of a grid square
outputpicture = "result.png"    # the name of the image to be created
outputtextfile = "result.txt"   # the name of a textfile, which contains the same info as the picture

VERSION = "0.1"
import os, sys, string
import Image        # external Image module - probably not provided in your python package

## my point class

class point3D:
	def __init__(self, number = "0", ordinate = 0, abszissa = 0, altitude = 0, code = "0"):
		self.n = number
		self.x = ordinate
		self.y = abszissa
		self.h = altitude
		self.c = code
	def __str__(self):
		return "%8s %6.3f %6.3f %3.3f %3s" % (self.n, self.x, self.y, self.h, self.c)


## minmax for pointlist

def minmax(pointlist):
	min_x = 10000000.0
	min_y = 10000000.0
	max_x = 0.0
	max_y = 0.0
	for punkt in pointlist:
		if punkt.x < min_x:
			min_x = punkt.x
		if punkt.x > max_x:
			max_x = punkt.x

		if punkt.y < min_y:
			min_y = punkt.y
		if punkt.y > max_y:
			max_y = punkt.y
	return min_x, min_y, max_x, max_y	


## determine grid

def griddef(gridspacing, min_x, min_y, max_x, max_y):
	grid_left =   min_x - (min_x % gridspacing)
	grid_bottom = min_y - (min_y % gridspacing)

	grid_top =   max_y + (gridspacing - (max_y % gridspacing))
	grid_right = max_x + (gridspacing - (max_x % gridspacing))

	rows = int((grid_top - grid_bottom) / gridspacing)
	columns = int((grid_right - grid_left) / gridspacing)

	return columns, rows, grid_left, grid_top, grid_right, grid_bottom 


## distribute points in grid

def distribute_points(punkteliste, gridspacing, columns, rows, grid_left, grid_top, grid_right, grid_bottom):
    grid = [[[]for x in range(rows)] for i in range(columns)]
     
    max_per_cell = 0
    
    while len(punkteliste) > 0:
    	point = punkteliste.pop()
    	row = int(((point.x - (point.x % gridspacing)) - grid_left) / gridspacing) 
    	column = int((grid_top - (point.y + (gridspacing - (point.y % gridspacing)))) / gridspacing)
    	grid[row][column].append(point)
    	
    	points_in_cell = len(grid[row][column]) 
        ######
    	if points_in_cell > max_per_cell:
    		max_per_cell = points_in_cell

    print "Maximum points in one cell: ", max_per_cell
    return grid, max_per_cell


## replace points in array field with their number

def count_points(grid):
    for column in range(columns):
    	for row in range(rows):
    		grid[column][row] = len(grid[column][row])
    return grid


## create picture from grid

def makepic(grid, outputpicture, max_per_cell):
	rows = len(grid[1])
	columns = len(grid)

	outputimage = Image.new('L', (columns, rows), (255))

	for row in range(columns):
		for column in range(rows):
			if grid[row][column] > 0:
				colour = 255 - int( grid[row][column] * 255 / max_per_cell) 
				outputimage.putpixel( (row,column), (colour) ) 

	outputimage.save(outputpicture, 'PNG')


## write data in file

def makefile(grid, textfilename):
    outputfile = file(textfilename, 'w')
    for row in range(columns):
    	outputfile.write("\n")
    	for column in range(rows):
    		outputfile.write(str(grid[row][column]) + " ")
    outputfile.close()

#############
# here basically "main" is starting
# read in the list of coordinates line by line, create coordinate objects
try:
	inputfile = file(sys.argv[1],'r')
except IOError:
	print "File ", sys.argv[1], " can not be opened for reading."

filecontent = inputfile.readlines()
inputfile.close()

pointlist = []

uncomplete_lines = 0
# FNORD
# This is the part, where the lines of the file are read in and parsed into pieces.
#
# The "if not ... " line will make the script ignore lines containing "!" or empty lines.
#
# The line "pieces = line.split()" splits the line divided by space(or tabs) in the array pieces[]. 
# You could put "pieces = line.split("-")", if the columns in your file are divided by "-".
#
# Most important is the line "pointlist.append(...". It creates a point3D object for each
# line by passing the the four attributes _number_, _easting_, _northing_ and _altitude_ to
# the constructor (+ a fifth attribute, the point-code, which is always set to "empty").
#
# Let us say your lines in the file look like this: "x-value, y-value, altitude, number"
# Then you would change the relevant line (as the values of pieces[] are set differently) to:
# pointlist.append(point3D(str(pieces[3]), float(pieces[0]), float(pieces[1]), float(pieces[2]), "empty",))
#
# If you ever did some programming yourself, you would not have needed this description.   :-)
# If you never programmed anyting yourself, it probably did not help you.   :-(

for line in filecontent:
	if not((line.find("!") == 0) or (line == "")):  # you might have to edit here
		try:
			pieces = line.split()                   # maybe here
			pointlist.append(point3D(str(pieces[0]), float(pieces[1]), float(pieces[2]), float(pieces[3]), "empty",)) # or here

		except IndexError:
			# print "Index out of range."
			uncomplete_lines += 1

output = str(len(pointlist)) + " points out of " 
output = output + sys.argv[1] + " read, with "
output = output + str(uncomplete_lines) + " uncomplete lines."
print output

# get border values
min_x, min_y, max_x, max_y = minmax(pointlist) 

print "smallest easting: ", min_x
print "largest easting:  ", max_x, " with delta: ", max_x - min_x
print "smallest northing:", min_y
print "largest northing: ", max_y, " with Delta: ", max_y - min_y

# determine gridparameters
columns, rows, grid_left, grid_top, grid_right, grid_bottom = griddef(gridspacing, min_x, min_y, max_x, max_y) 

print "lines: %s  rows: %s with %s m grid --> %s cells" % (rows, columns, gridspacing, rows*columns)

# points in list are being distributed over grid
grid, max_per_cell = distribute_points(pointlist, gridspacing, columns, rows, grid_left, grid_top, grid_right, grid_bottom)

# change points in quantitiy
grid = count_points(grid)

# output as textfile
makefile(grid, outputtextfile);

# output as picture
makepic(grid, outputpicture, max_per_cell)