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
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.