Search code examples
pythonarraysmultidimensional-arrayshapely

Apply function to 2D array with x,y coordinates as parameters


I am planning to highlight an area in elevation data in 3D as explained [here][1]. Therefore, I need to have an array where all the points are set to 0 that should not be highlighted. The data array looks like below:

y          51.380591   51.380769  ...   51.408314   51.408493
x                                 ...                        
9.171581  223.699997  228.199997  ...  297.500000  297.899994
9.171868  220.199997  223.199997  ...  297.600006  297.600006
9.172156  218.900009  220.600006  ...  296.100006  296.100006
9.172443  218.000000  218.100006  ...  295.000000  295.000000
9.172731  223.400009  220.699997  ...  297.100006  297.100006
             ...         ...  ...         ...         ...
9.212988  252.600006  254.100006  ...  313.100006  313.000000
9.213276  257.399994  260.700012  ...  313.600006  313.200012
9.213563  259.500000  262.700012  ...  314.700012  314.000000
9.213851  261.399994  264.600006  ...  315.800018  315.399994
9.214139  262.500000  265.800018  ...  317.399994  317.300018

[149 rows x 157 columns]

My goal is to only keep values that fulfill a condition based on their x,y location. All other values should be set to 0. In detail I would like to define a 2D polygon. Every point that is not within this polygon should be set to 0.

An example polygon looks like this.

coords = [(9.185, 51.38), (9.175, 51.385), (9.175, 51.4), (9.2, 51.395)]
poly = Polygon(coords)

The function to check if a point is within a polygon is provided by [shapely][2].

from shapely.geometry import Point, Polygon
Point(x,y).within(poly)

So the function that should be applied to all cells looks like this:

def KeepValue(x,y,z,poly):
    if (Point(x,y).within(poly)):
        return z
    return 0

My question is Which function should I use to apply the within function to all cells? The major problem for me is to find a way to hand over especially the x,y and z parameters.

I have tried apply/applymap but without good results. [1]: How to mark an area in plotly 3D surface plot? [2]: https://automating-gis-processes.github.io/CSC18/lessons/L4/point-in-polygon.html


Solution

  • I have not found a good solution so far, but iterating the columns and then the values is working properly:

    contour_data = pd.read_csv(r"C:\Elevation.xyz", delimiter=' ', names=["x","y","z"])
    Z = contour_data.pivot_table(index='x', columns='y', values='z').T
    X_unique = np.sort(contour_data.x.unique())
    Y_unique = np.sort(contour_data.y.unique())
    
    coords = [(9.185, 51.39), (9.175, 51.39), (9.175, 51.4), (9.2, 51.395)]
    poly = Polygon(coords)
    
    marked_area=Z.copy()
    i=0
    for x in X_unique:
        j=0
        for z in Z.iloc[i]:
            if (Point(x,Y_unique[j]).within(poly)):
                marked_area.iloc[i,j]=z+0.1
            else:
                marked_area.iloc[i,j]=0
            j=j+1
        i=i+1
    

    This may not be the most efficient approach, but this fits my needs.