3
\$\begingroup\$

I have a set of data which represents a set of nodes, each node is associated with a value (represented by a color in the image). What I want to achieve is selectively changing those values.

The mesh represents a porous system (say a rock for example) model. The pressure in my system is specified at the nodes. I want to be able to assign initial conditions for the pressure at specific nodes only (the ones located inside the polygon). So the weights of my nodes are the pressure at that node.

I actually want to define a polygon, and attribute a value to each vertex (think of it as a weight), and using the weight of the vertex and the distance from the vertices to each node INSIDE the polygon to correct the value for that node.

This is what my output look like:

enter image description here

I am working on an algorithm which takes in a set of values of the form [x,y,z] and another with [value,value,value,value]. Both have the same number of rows. e.i, the rows in the first input are the location of the node, and the rows in the second the values associated with that node.

I made an algorithm which takes in a set of points forming a polygon and a set of weights corresponding to each vertex of that polygon.

I then scan my merged inputs and replace the value of any node found to be inside the polygon. The value is defined by the algorithm in this post. If the node is not inside the polygon it's value is preserved.

I then write all the new_values to a file.

My concern is that I will not be able to handle large inputs of a few million nodes without a MemoryError. At the moment handling a 9239 rows of inputs, takes me 9 seconds.

This is my code:

# for PIP problem
import shapely.geometry as shapely
# for plot
import matplotlib.pyplot as plt
# for handling data
import csv
import itertools
# for timing
import time

#=================================================================================
                # POINT IN POLYGONE PROBLEM
#=================================================================================

class MyPoly(shapely.Polygon):
    def __init__(self,points):
        closed_path = list(points)+[points[0]]
        super(MyPoly,self).__init__(closed_path)
        self.points = closed_path
        self.points_shapely = [shapely.Point(p[0],p[1]) for p in closed_path]

def convert_to_shapely_points_and_poly(poly,points):
    poly_shapely = MyPoly(poly)
    points_shapely = (shapely.Point(p[0],p[1]) for p in points)
    return poly_shapely,points_shapely    


def isBetween(a, b, c): #is c between a and b ?
    crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
    if abs(crossproduct) > 0.01 : return False   # (or != 0 if using integers)

    dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y)
    if dotproduct < 0 : return False

    squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
    if dotproduct > squaredlengthba: return False

    return True

def get_edges(poly):
    # get edges
    edges = []
    for i in range(len(poly.points)-1):
        t = [poly.points_shapely[i],poly.points_shapely[i+1]]
        edges.append(t)
    return edges


def inPoly(poly,point, inclusive):
    if poly.contains(point) == True:
        return 1
    elif inclusive:
        for e in get_edges(poly):
            if isBetween(e[0],e[1],point):
                return 1
    return 0

def plot(poly_init,points_init, inclusive = True):
    #convert to shapely poly and points
    poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)

    #plot polygon
    plt.plot(*zip(*poly.points))

    #plot points
    xs,ys,cs = [],[],[]
    for point in points:
        xs.append(point.x)
        ys.append(point.y)
        color = inPoly(poly,point, inclusive)
        cs.append(color)
        print point,":", color
    plt.scatter(xs,ys, c = cs , s = 20*4*2)

    #setting limits
    axes = plt.gca()
    axes.set_xlim([min(xs)-5,max(xs)+50])
    axes.set_ylim([min(ys)-5,max(ys)+10])

    plt.show()


# TESTS ========================================================================
#set up poly
polys = {
1 : [[10,10],[10,50],[50,50],[50,80],[100,80],[100,10]], # test rectangulary shape
2 : [[20,10],[10,20],[30,20]], # test triangle
3 : [[0,0],[0,10],[20,0],[20,10]], # test bow-tie
4 : [[0,0],[0,10],[20,10],[20,0]], # test rect clockwise
5 : [[0,0],[20,0],[20,10],[0,10]] # test rect counter-clockwise
}

#points to check
points = {
1 : [(10,25),(50,75),(60,10),(20,20),(20,60),(40,50)], # rectangulary shape test pts
2 : [[20,10],[10,20],[30,20],[-5,0],[20,15]] , # triangle  test pts
3 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,5]],  # bow-tie shape test pts
4 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]],  # rect shape test pts
5 : [[0,0],[0,10],[20,0],[20,10],[10,0],[10,5],[15,2],[30,8]]  # rect shape test pts
}

for data in zip(polys.itervalues(),points.itervalues()):
    plot(data[0],data[1], True)

#================================================================================
#                           WEIGHTING FUNCTION
#================================================================================

def add_weights(poly, weights):
    poly.weights = [float(w) for w in weights]+[weights[0]] #need to add the first weight
                                                          # at the end to account for
                                                          # the first point being added to close the loop

def distance(a,b):
    dist = ( (b.x - a.x)**2 + (b.y - a.y)**2 )**0.5
    if dist == 0: dist = 0.000000001 
    return dist

def get_weighted_sum(poly, point):
    return sum([poly.weights[n]/distance(point,p) for n,p in enumerate(poly.points_shapely) if poly.weights[n] != 'nan'])

def get_weighted_dist(poly, point):
    return sum([1/distance(point,p) for n,p in enumerate(poly.points_shapely) if poly.weights[n] != 'nan'])

def get_point_weighted_value(poly, point):
    return get_weighted_sum(poly,point)/get_weighted_dist(poly,point)

#==============================================================================
    #           GETTING THE DATA inside the Polygone
#==============================================================================

'''Function Definitions'''

def data_extraction(filename,start_line,node_num,span_start,span_end):
    with open(filename, "r") as myfile:
        file_= csv.reader(myfile, delimiter=' ')  #extracts data from .txt as lines
        return (x for x in [filter(lambda a: a != '', row[span_start:span_end]) \
        for row in itertools.islice(file_, start_line, node_num)])

def merge_data(msh_data,reload_data):
    return (zip(msh_data,reload_data))


def edit_value(data, poly_init, weights):
    #make x,y coordinates of the data into points
    points_init = ([float(pair[0][0]),float(pair[0][2])]for pair in data)
    #convert to shapely poly and points
    poly,points = convert_to_shapely_points_and_poly(poly_init,points_init)
    add_weights(poly,weights)
    #fliter out points in polygon
    new_pair = []
    for n,point in enumerate(points):
        if inPoly(poly, point, True):
            value = str(get_point_weighted_value(poly, point))
            new_pair.append([value,value,value,value])
        else:
            new_pair.append(data[n][1])

    return (x for x in new_pair)


def make_file(filename,data):
    with open(filename, "w") as f:
        f.writelines('\t'.join(i) + '\n' for i in data)
        f.close()

def run(directory_path,poly_list,weight_list):
    ''' Directory and File Names'''
    msh_file = '\\nodes.txt'
    reload_file = '\\values.txt'
    new_reload_file = '\\new_values.txt'
    dir_path = directory_path
    msh_path = dir_path+msh_file
    reload_path = dir_path+reload_file

    '''Running the Code with your data'''
    mesh_data = data_extraction(msh_path,0,54,1,4)
    reload_data = data_extraction(reload_path,0,54,0,7)
    data = merge_data(mesh_data,reload_data)
    new_values = edit_value(data,poly_list,weight_list)
    make_file(dir_path+new_reload_file,new_values)        



t0 = time.time()
run("M:\\MyDocuments"
    ,([75,-800],[50,-900],[50,-1350],[90,-1000],[100,-900])
    ,(5.0e5,1e6,1e8,5.0e5,1.0e7))
t1= time.time()
print t1-t0

And this is sample data for you to play with (you will need to save it as a values.txt in the folder of interest):

1.067896706746556e+006  8.368595971460000e+006  1.068728658407457e+006  8.368595971460000e+006  
2.844581459224940e+005  8.613334125294963e+006  2.846631849296530e+005  8.613337865004616e+006  
1.068266636556349e+006  8.368595971460000e+006  1.069097067800019e+006  8.368595971460000e+006  
2.844306728256134e+005  8.613334269264592e+006  2.846366503960088e+005  8.613338015263893e+006  
2.646871122251647e+003  9.280390372578276e+006  2.647124079593603e+003  9.279361848151272e+006  
2.645513962411728e+003  9.280388336827532e+006  2.645732877622660e+003  9.279359747270351e+006  
1.067996132697676e+006  8.368595971460000e+006  1.068827019510901e+006  8.368595971460000e+006  
1.068040363056876e+006  8.368595971460000e+006  1.068870797759632e+006  8.368595971460000e+006  
1.068068562573701e+006  8.368595971460000e+006  1.068898735336173e+006  8.368595971460000e+006  
1.068088894288788e+006  8.368595971460000e+006  1.068918905897983e+006  8.368595971460000e+006  
1.068104407561180e+006  8.368595971460000e+006  1.068934323713974e+006  8.368595971460000e+006  
1.068116587634527e+006  8.368595971460000e+006  1.068946455001944e+006  8.368595971460000e+006  
1.068126287610437e+006  8.368595971460000e+006  1.068956140100951e+006  8.368595971460000e+006  
1.068134059350058e+006  8.368595971460000e+006  1.068963921144020e+006  8.368595971460000e+006  
1.068140293705664e+006  8.368595971460000e+006  1.068970181121935e+006  8.368595971460000e+006  
1.068145286907994e+006  8.368595971460000e+006  1.068975209852979e+006  8.368595971460000e+006  
1.068149274285654e+006  8.368595971460000e+006  1.068979237556288e+006  8.368595971460000e+006  
1.068152448234754e+006  8.368595971460000e+006  1.068982452745734e+006  8.368595971460000e+006  
1.068154968237062e+006  8.368595971460000e+006  1.068985012157432e+006  8.368595971460000e+006  
1.068156966832556e+006  8.368595971460000e+006  1.068987046589365e+006  8.368595971460000e+006  
1.068158553584154e+006  8.368595971460000e+006  1.068988664694503e+006  8.368595971460000e+006  
1.068159818109515e+006  8.368595971460000e+006  1.068989955820237e+006  8.368595971460000e+006  
1.068160832713088e+006  8.368595971460000e+006  1.068990992449035e+006  8.368595971460000e+006  
1.068161654841110e+006  8.368595971460000e+006  1.068991832481593e+006  8.368595971460000e+006  
1.068162329417389e+006  8.368595971460000e+006  1.068992521432464e+006  8.368595971460000e+006  
1.068162891039026e+006  8.368595971460000e+006  1.068993094522149e+006  8.368595971460000e+006  
1.068163365980015e+006  8.368595971460000e+006  1.068993578612505e+006  8.368595971460000e+006  
1.068163773959426e+006  8.368595971460000e+006  1.068993993936813e+006  8.368595971460000e+006  
1.068164129647837e+006  8.368595971460000e+006  1.068994355591033e+006  8.368595971460000e+006  
1.068164443906155e+006  8.368595971460000e+006  1.068994674772899e+006  8.368595971460000e+006  
1.068164724771358e+006  8.368595971460000e+006  1.068994959776965e+006  8.368595971460000e+006  
1.068164978220522e+006  8.368595971460000e+006  1.068995216771951e+006  8.368595971460000e+006  
1.068165208742598e+006  8.368595971460000e+006  1.068995450386572e+006  8.368595971460000e+006  
1.068165419759155e+006  8.368595971460000e+006  1.068995664143054e+006  8.368595971460000e+006  
1.068165613928702e+006  8.368595971460000e+006  1.068995860772074e+006  8.368595971460000e+006  
1.068165793358832e+006  8.368595971460000e+006  1.068996042433189e+006  8.368595971460000e+006  
1.068165959755557e+006  8.368595971460000e+006  1.068996210870328e+006  8.368595971460000e+006  
1.068166114528858e+006  8.368595971460000e+006  1.068996367521690e+006  8.368595971460000e+006  
1.068166258863683e+006  8.368595971460000e+006  1.068996513593785e+006  8.368595971460000e+006  
1.068166393771001e+006  8.368595971460000e+006  1.068996650114520e+006  8.368595971460000e+006  
1.068166520124043e+006  8.368595971460000e+006  1.068996777970799e+006  8.368595971460000e+006  
1.068166638684427e+006  8.368595971460000e+006  1.068996897935547e+006  8.368595971460000e+006  
1.068166750121474e+006  8.368595971460000e+006  1.068997010687638e+006  8.368595971460000e+006  
1.068166855027363e+006  8.368595971460000e+006  1.068997116827437e+006  8.368595971460000e+006  
1.068166953929060e+006  8.368595971460000e+006  1.068997216889020e+006  8.368595971460000e+006  
1.068167047297063e+006  8.368595971460000e+006  1.068997311349124e+006  8.368595971460000e+006  
1.068167135553131e+006  8.368595971460000e+006  1.068997400635036e+006  8.368595971460000e+006  
1.068167219077452e+006  8.368595971460000e+006  1.068997485131862e+006  8.368595971460000e+006  
1.068167298214444e+006  8.368595971460000e+006  1.068997565188462e+006  8.368595971460000e+006  
1.068167373276848e+006  8.368595971460000e+006  1.068997641121696e+006  8.368595971460000e+006  
1.068167444552389e+006  8.368595971460000e+006  1.068997713223352e+006  8.368595971460000e+006  
1.068167512315698e+006  8.368595971460000e+006  1.068997781772706e+006  8.368595971460000e+006  
1.068167576851623e+006  8.368595971460000e+006  1.068997847061655e+006  8.368595971460000e+006  
1.068167638524622e+006  8.368595971460000e+006  1.068997909469268e+006  8.368595971460000e+006  
1.068167697990453e+006  8.368595971460000e+006  1.068997969688582e+006  8.368595971460000e+006 

and a sample of the node data (will need to be saved as nodes.txt):

0 0.26 0 -800.0
1 0.26 0 -1062.5
2 143.0 0 -800.0
3 143.0 0 -1062.5
4 0.26 0 -1150.0
5 143.0 0 -1150.0
6 1.17057404659 0 -800.0
7 2.10837283486 0 -800.0
8 3.07421037484 0 -800.0
9 4.06892500937 0 -800.0
10 5.0933801161 0 -800.0
11 6.14846485319 0 -800.0
12 7.23509501358 0 -800.0
13 8.35421377171 0 -800.0
14 9.50679247207 0 -800.0
15 10.6938315019 0 -800.0
16 11.9163611811 0 -800.0
17 13.1754426445 0 -800.0
18 14.472168711 0 -800.0
19 15.8076649025 0 -800.0
20 17.1830903981 0 -800.0
21 18.59963901 0 -800.0
22 20.0585402542 0 -800.0
23 21.5610604197 0 -800.0
24 23.1085036396 0 -800.0
25 24.7022130583 0 -800.0
26 26.3435719675 0 -800.0
27 28.0340049986 0 -800.0
28 29.7749793906 0 -800.0
29 31.5680062618 0 -800.0
30 33.4146418792 0 -800.0
31 35.3164890883 0 -800.0
32 37.2751986287 0 -800.0
33 39.2924705661 0 -800.0
34 41.3700558588 0 -800.0
35 43.5097577902 0 -800.0
36 45.7134335243 0 -800.0
37 47.9829957969 0 -800.0
38 50.3204145133 0 -800.0
39 52.7277184748 0 -800.0
40 55.2069971476 0 -800.0
41 57.7604024715 0 -800.0
42 60.3901507176 0 -800.0
43 63.0985244223 0 -800.0
44 65.8878743782 0 -800.0
45 68.7606216458 0 -800.0
46 71.7192596577 0 -800.0
47 74.7663564153 0 -800.0
48 77.904556725 0 -800.0
49 81.1365844387 0 -800.0
50 84.4652448319 0 -800.0
51 87.8934270922 0 -800.0
52 91.4241067682 0 -800.0
53 95.0603483625 0 -800.0
54 98.8053080549 0 -800.0
55 102.662236327 0 -800.0
\$\endgroup\$
  • \$\begingroup\$ Could you come up with a less generic title that actually explains what this code does? Also, you've told us a lot about how you're manipulating the data, but you've explained nothing about what you are really trying to accomplish. \$\endgroup\$ – 200_success Sep 16 '16 at 16:51
  • \$\begingroup\$ I have edited the question with a visual representation of what I want to achieve. I also modified the title. \$\endgroup\$ – Sorade Sep 16 '16 at 17:02
  • \$\begingroup\$ @Sorade: Can you say more about the motivation for this problem? What do these nodes represent? What does the polygon represent? Why do you want to change the nodes inside the polygon? What is the significance of the weighted distance to the vertices of the polygon? \$\endgroup\$ – Gareth Rees Sep 18 '16 at 9:10
  • 1
    \$\begingroup\$ I've just added some information in the question. The mesh is a model of a porous system (say a rock) and the node weight is the pressure at the node. \$\endgroup\$ – Sorade Sep 18 '16 at 9:22
2
\$\begingroup\$

I find that I get a 7× speedup simply by rewriting the inPoly function as follows:

def inPoly(poly, point, inclusive=False):
    """Return true if point is inside the polygon, or inclusive is true
    and point is on the boundary of the polygon.

    """
    return poly.contains(point) or inclusive and poly.boundary.contains(point)

See the shapely documentation for BaseGeometry.boundary for an explanation of how this works.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.