import collections
import json
import math
from multiprocessing import Process, Queue
import numpy as np
from scipy.spatial import distance
from shapely.geometry import MultiPolygon, Point, Polygon
[docs]def worker(input, output):
for func, args in iter(input.get, "STOP"):
result = func(*args)
output.put(result)
[docs]def maxDist1(p, scale):
coords = np.array(p) * scale
array = distance.cdist(coords, coords, "euclidean")
shape = array.shape
max = []
for i in range(shape[0]):
max.append(np.max(array[i]))
max_of_max = np.max(max)
return max_of_max
[docs]def simple_distance(p1, p2, scale):
x0 = (p1[0] - p2[0]) * scale
y0 = math.fabs((p1[1] - p2[1]) * scale)
y0 = (p1[1] - p2[1]) * scale
tan = math.degrees(math.atan2(x0, y0))
dist = math.sqrt(x0 * x0 + y0 * y0)
return x0 * x0 + y0 * y0, tan
# Function to find the maximum
# distance between any two points
[docs]def maxDist(p, scale):
n = len(p)
# print ('length of p ', n)
maxm = 0
tan = 0
p_0 = p[0]
p_1 = p[1]
# Iterate over all possible pairs
for i in range(n):
for j in range(i + 1, n):
# Update maxm
dist, tan1 = simple_distance(p[i], p[j], scale)
if dist > maxm:
maxm = dist
tan = tan1
p_0 = p[i]
p_1 = p[j]
tan = 90.0 + tan
if tan < 0.0:
tan = tan + 180.0
if tan > 180.0:
tan = tan - 180.0
return (math.sqrt(maxm), tan, p_0, p_1)
[docs]def process_json_file(in_data, pixel_size=0.0):
num_contours = in_data["num_contours"]
# print('num_contours', num_contours)
try:
coords = in_data["coords"]
except:
coords = []
if len(coords) == 0 and in_data["manual"] == False:
polygon_list = []
return polygon_list, coords
# print('number of separate contours', num_contours)
outer_list = []
for j in range(num_contours):
poly_coord = in_data[str(j)]
if len(poly_coord) > 2:
p = Polygon(poly_coord)
inner_list = [j, p.area]
outer_list.append(inner_list)
length = len(outer_list)
if length > 0:
arr2d = np.array(outer_list)
columnIndex = 1
sortedArr = arr2d[arr2d[:, columnIndex].argsort()[::-1]] # sorts in ascending order
# print('sorted arr2d', sortedArr)
# find n largest values
# print('**** coords', coords)
if in_data["manual"]:
n = 1
else:
n = len(coords)
polygon_list = []
las = []
las_pa = []
p_max = []
for i in range(n):
# print('outer iteration', i)
found = False
# print('coordinate ', coords[i])
if not in_data["manual"]:
Pt = Point(coords[i])
else:
Pt = Point([0.0, 0.0])
for l in range(length):
# print('inner iteration', l)
rslt = sortedArr[l]
contour_number = int(rslt[0])
poly_coord = in_data[str(contour_number)]
p = Polygon(poly_coord)
if len(p.exterior.coords) > 500:
old = len(p.exterior.coords)
p = p.simplify(0.5, True)
new = len(p.exterior.coords)
if new < old:
in_data[str(contour_number)] = list(p.exterior.coords)
if p.contains(Pt) or in_data["manual"]:
# print('********* containing contour', l)
polygon_list.append(p)
if pixel_size > 0.0:
# result = maxDist1(poly_coord,pixel_size)
# print('maxDist1 gives ', result)
# las.append(result)
result = maxDist(poly_coord, pixel_size)
las.append(result[0])
las_pa.append(result[1])
p_max.append(result[2])
p_max.append(result[3])
if pixel_size > 0:
result = maxDist(p_max, pixel_size)
max_las = result[0]
max_pa = result[1]
# print('Source has maximum angular size and position angle', max_las, max_pa)
return polygon_list, coords, max_las, las, max_pa, las_pa
else:
return polygon_list, coords
[docs]def get_contained_points(y, xmin, xmax, i, p):
x_list = []
y_list = []
for x in range(xmin, xmax):
pt = Point(x, y)
# Testing for contained points can be somewhat slow if you have
# a big polygon containing many points
# Note that we flip x and y coordinates going from matplotlib
# display coordinates to numpy array coordinates
if p.contains(pt):
x_list.append(y)
y_list.append(x)
if len(x_list) != len(y_list):
# print('error x_list and y_list have different lengths', len(x_list), len(y_list))
return (-1, -1, -1.0, -1.0)
else:
if len(x_list) > 0 and len(y_list) > 0:
return (i, y, (x_list, y_list))
else:
return (-1, -1, -1.0, -1.0)
[docs]def get_interior_locations(polygon_list):
# print('Getting points interior to polygon boundaries ...')
num_processors = 1
# if num_processors <= 2 and len(polygon_list) > 1:
if num_processors <= 2:
try:
import multiprocessing
processors = multiprocessing.cpu_count()
if processors > num_processors:
num_processors = processors
except:
pass
# print ('*** setting final number of processors to',num_processors)
TASKS = []
# print('getting pixels inside polygon, working ...')
for i in range(len(polygon_list)):
p = polygon_list[i]
# Can we simplify the polygon?
if len(p.exterior.coords) > 500:
old = len(p.exterior.coords)
p = p.simplify(0.5, True)
new = len(p.exterior.coords)
if new < old:
# print('simplify shrank polygon points from ', old, ' to ', new)
polygon_list[i] = p
(minx, miny, maxx, maxy) = p.bounds
# set up parallelprocessing
for y in range(int(miny - 1), int(maxy + 1)):
TASKS.append((get_contained_points, (y, int(minx - 1), int(maxx + 1), i, p)))
NUMBER_OF_PROCESSES = num_processors
# Create queues
task_queue = Queue()
done_queue = Queue()
# Submit tasks
for task in TASKS:
task_queue.put(task)
# Start worker processes
for i in range(NUMBER_OF_PROCESSES):
Process(target=worker, args=(task_queue, done_queue)).start()
interior_polygons = {}
result_sum = 0
first = True
for i in range(len(polygon_list)):
interior_polygons[i] = ([], [])
# get output from tasks
for i in range(len(TASKS)):
result = done_queue.get(timeout=2000)
if result[0] >= 0:
data = result[2]
if len(data[0]) != len(data[1]):
print("retrieved data have different lengths!")
else:
data = interior_polygons[result[0]]
x_list = data[0]
y_list = data[1]
x_list = x_list + result[2][0]
y_list = y_list + result[2][1]
interior_polygons[result[0]] = (x_list, y_list)
# make sure things shut down
for i in range(num_processors):
task_queue.put("STOP")
interior_points = []
for i in range(len(polygon_list)):
data = interior_polygons[i]
x_list = data[0]
y_list = data[1]
interior_points.append((i, (np.array(x_list), np.array(y_list))))
return interior_points