#! /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)