'Triangle interpolation script not working
I solved a part of my question! So I removed the original question and updated it to this one:
I am making a wind interpolation script. There are several weather stations that surround the Wadden Sea. I have done measurements at a certain location (P). And I have a .csv file with all the wind speeds, per station, per time, that looks like this:
STN DT HH FF
0 229 2-1-2021 9 20.0
1 229 2-1-2021 10 20.0
2 229 2-1-2021 11 10.0
3 229 2-1-2021 12 10.0
4 229 2-1-2021 13 30.0
... ... .. ...
62747 277 12-3-2022 20 50.0
62748 277 12-3-2022 21 50.0
62749 277 12-3-2022 22 50.0
62750 277 12-3-2022 23 60.0
62751 277 12-3-2022 24 50.0
To derive the wind speed at point (P), I will use a triangle interpolation. The script try to fit (P) within a triangle that uses the possible stations as vertices. This works fine now, after some help. I got it in this script:
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 16 09:59:31 2022
@author: ###
"""
import itertools
import math
from pyproj import Transformer
import numpy as np
from math import sqrt
import os
import pandas as pd
from pathlib import Path
from matplotlib import pyplot as plt
## data
os.chdir('C:\\Users\\###\\sensoren\\data')
os.getcwd()
wind = pd.read_csv('data_totaal_2.csv', usecols= [ 0, 1,2 , 5],skiprows = 33, names = ['STN','DT','HH','FF'])
## closest point
def min_distance(x, y, iterable):
list_of_distances = list(map(lambda t: sqrt(pow(t[0]-x,2)+pow(t[1]-y,2)),iterable))
min_res = min(list_of_distances)
index_of_min = list_of_distances.index(min_res)
return iterable[index_of_min]
## check if point is inside triangle
def area(x1, y1, x2, y2, x3, y3):
return abs((x1 * (y2 - y3) +
x2 * (y3 - y1) +
x3 * (y1 - y2)) / 2.0)
def is_inside(x1, y1, x2, y2, x3, y3, x, y):
A = area(x1, y1, x2, y2, x3, y3)
A1 = area(x, y, x2, y2, x3, y3)
A2 = area(x1, y1, x, y, x3, y3)
A3 = area(x1, y1, x2, y2, x, y)
return A == A1 + A2 + A3
## calculate perimeter
def perimeter(points):
x1, y1, x2, y2, x3, y3 = *points[0], *points[1], *points[2]
return (math.hypot(x1-x2, y1-y2) +
math.hypot(x2-x3, y2-y3) +
math.hypot(x3-x1, y3-y1))
## barycentric interpolation
def bary_interpol (inter_cor_x1, inter_cor_x2, inter_cor_x3, inter_cor_y1,
inter_cor_y2, inter_cor_y3, inter_cor_px, inteR_cor_py,
inter_wind_v1, inter_wind_v2, inter_wind_v3):
W_v1 = ((inter_cor_y2-inter_cor_y3)*(inter_cor_px-inter_cor_x3)+
(inter_cor_x3-inter_cor_x2)*(inter_cor_py-inter_cor_y3))/((inter_cor_y2-inter_cor_y3)*(inter_cor_x1-inter_cor_x3)+
(inter_cor_x3-inter_cor_x2)*(inter_cor_y1-inter_cor_y3))
W_v2 = ((inter_cor_y3-inter_cor_y1)*(inter_cor_px-inter_cor_x3)+
(inter_cor_x1-inter_cor_x3)*(inter_cor_py-inter_cor_y3))/((inter_cor_y2-inter_cor_y3)*(inter_cor_x1-inter_cor_x3)+
(inter_cor_x3-inter_cor_x2)*(inter_cor_y1-inter_cor_y3))
W_v3 = 1 - W_v1 - W_v2
inter_wind_pv = (W_v1 * inter_wind_v1) + (W_v2 * inter_wind_v2) + (W_v3 * inter_wind_v3)
return inter_wind_pv
## datum shift
transformer = Transformer.from_crs("epsg:4326", "epsg:28992")
RD_229 = transformer.transform(53.04309, 4.81969) #texelhors_229
RD_235 = transformer.transform(52.91382, 4.79498) #de_kooy_235
RD_242 = transformer.transform(53.25337, 4.94041) #vlieland_242
RD_251 = transformer.transform(53.39335, 5.34623) #terschelling_251
RD_267 = transformer.transform(52.89709, 5.38331) #stavoren_267
RD_277 = transformer.transform(53.40848, 6.19614) #lauwersoog_277
P = transformer.transform(53.02, 5.02) #location 1
# =============================================================================
# P = transformer.transform(50.00, 5.00) #location outside any triangle
# =============================================================================
print(RD_229)
print(RD_235)
print(RD_242)
print(RD_251)
print(RD_267)
print(RD_277)
STN_229_x = RD_229[0]
STN_235_x = RD_235[0]
STN_242_x = RD_242[0]
STN_251_x = RD_251[0]
STN_267_x = RD_267[0]
STN_277_x = RD_277[0]
points = [(RD_229[0],RD_229[1]),(RD_235[0],RD_235[1]),(RD_242[0],RD_242[1]),(RD_251[0],RD_251[1]),(RD_267[0],RD_267[1]),(RD_277[0],RD_277[1])]
## find best triangle or point
triangle_fit = []
for triangle in itertools.combinations(points, 3):
p1, p2, p3 = triangle
if is_inside(*p1, *p2, *p3, *P):
triangle_fit.append(triangle)
triangle_fit.sort(key=perimeter)
# Sort by perimeter.
triangle_fit.append(triangle)
triangle_fit.sort(key=perimeter)
inter_cor = triangle_fit[0]
inter_cor_x1 = inter_cor[0][0]
inter_cor_x2 = inter_cor[1][0]
inter_cor_x3 = inter_cor[2][0]
inter_cor_y1 = inter_cor[0][1]
inter_cor_y2 = inter_cor[1][1]
inter_cor_y3 = inter_cor[2][1]
inter_cor_px = P[0]
inter_cor_py = P[1]
print(inter_cor_x1)
print(inter_cor_x2)
print(inter_cor_x3)
print(inter_cor_y1)
print(inter_cor_y2)
print(inter_cor_y3)
if inter_cor_x1 == STN_229_x:
STN_v1 = 229
elif inter_cor_x1 == STN_235_x:
STN_v1 = 235
elif inter_cor_x1 == STN_242_x:
STN_v1 = 242
elif inter_cor_x1 == STN_251_x:
STN_v1 = 251
elif inter_cor_x1 == STN_267_x:
STN_v1 = 267
elif inter_cor_x1 == STN_277_x:
STN_v1 = 277
else:
print("no matching stations in v1")
if inter_cor_x2 == STN_229_x:
STN_v2 = 229
elif inter_cor_x2 == STN_235_x:
STN_v2 = 235
elif inter_cor_x2 == STN_242_x:
STN_v2 = 242
elif inter_cor_x2 == STN_251_x:
STN_v2 = 251
elif inter_cor_x2 == STN_267_x:
STN_v2 = 267
elif inter_cor_x2 == STN_277_x:
STN_v2 = 277
else:
print("no matching stations in v2")
if inter_cor_x3 == STN_229_x:
STN_v3 = 229
elif inter_cor_x3 == STN_235_x:
STN_v3 = 235
elif inter_cor_x3 == STN_242_x:
STN_v3 = 242
elif inter_cor_x3 == STN_251_x:
STN_v3 = 251
elif inter_cor_x3 == STN_267_x:
STN_v3 = 267
elif inter_cor_x3 == STN_277_x:
STN_v3 = 277
else:
print("no matching stations in v1")
inter_wind_v1 = wind[wind['STN'] == STN_v1]
inter_wind_v2 = wind[wind['STN'] == STN_v2]
inter_wind_v3 = wind[wind['STN'] == STN_v3]
df_add = inter_wind_v1.FF.add(inter_wind_v1.FF, fill_value=0)
inter_wind_v1 = wind[wind['STN'] == STN_v1].reset_index()
# =============================================================================
# print(f"inter_wind_v1:\n{inter_wind_v1}")
# =============================================================================
inter_wind_v2 = wind[wind['STN'] == STN_v2].reset_index()
# =============================================================================
# print(f"inter_wind_v2:\n{inter_wind_v2}")
# =============================================================================
inter_wind_v3 = wind[wind['STN'] == STN_v3].reset_index()
# =============================================================================
# print(f"inter_wind_v3:\n{inter_wind_v3}")
# =============================================================================
interp = pd.DataFrame(inter_wind_v1['FF']).rename(columns = {'FF':'inter_wind_v1'})
interp['inter_wind_v2'] = inter_wind_v2['FF']
interp['inter_wind_v3'] = inter_wind_v3['FF']
# =============================================================================
# print(interp)
# =============================================================================
def doInterp(inter_wind_v1, inter_wind_v2, inter_wind_v3):
return bary_interpol (inter_cor_x1, inter_cor_x2, inter_cor_x3, inter_cor_y1,
inter_cor_y2, inter_cor_y3, inter_cor_px, inter_cor_py,
inter_wind_v1, inter_wind_v2, inter_wind_v3)
interp['inter_wind_pv'] = interp.apply(lambda x: doInterp(x['inter_wind_v1'], x['inter_wind_v2'], x['inter_wind_v3']), axis=1)
print(interp)
filepath = Path('C:/Users/###/sensoren/data/out.csv')
filepath.parent.mkdir(parents=True, exist_ok=True)
interp.to_csv(filepath)
os.makedirs('C:/Users/###/sensoren/data/', exist_ok=True)
interp.to_csv('C:/Users/###/sensoren/data/out2.csv')
data_s = np.array([
[RD_229],
[RD_235],
[RD_242],
[RD_251],
[RD_267],
[RD_277]
])
data_p = np.array([
[P]
])
data_i = np.array([
[ inter_cor_x1, inter_cor_y1],
[ inter_cor_x2, inter_cor_y2],
[ inter_cor_x3, inter_cor_y3]
])
x_s, y_s = data_s.T
x_p, y_p = data_p. T
x_i, y_i = data_i. T
plt.scatter( x_p, y_p, c="red")
plt.scatter( x_s, y_s, c="blue")
plt.scatter( x_i, y_i, c="green")
plt.show()
## COPY FOR POINT 2 ##
P = transformer.transform(53.29, 5.40) #location 2
## find best triangle or point
triangle_fit = []
for triangle in itertools.combinations(points, 3):
p1, p2, p3 = triangle
if is_inside(*p1, *p2, *p3, *P):
triangle_fit.append(triangle)
triangle_fit.sort(key=perimeter)
#sort by perimeter
triangle_fit.append(triangle)
triangle_fit.sort(key=perimeter)
inter_cor = triangle_fit[0]
inter_cor_x1 = inter_cor[0][0]
inter_cor_x2 = inter_cor[1][0]
inter_cor_x3 = inter_cor[2][0]
inter_cor_y1 = inter_cor[0][1]
inter_cor_y2 = inter_cor[1][1]
inter_cor_y3 = inter_cor[2][1]
inter_cor_px = P[0]
inter_cor_py = P[1]
if inter_cor_x1 == STN_229_x:
STN_v1 = 229
elif inter_cor_x1 == STN_235_x:
STN_v1 = 235
elif inter_cor_x1 == STN_242_x:
STN_v1 = 242
elif inter_cor_x1 == STN_251_x:
STN_v1 = 251
elif inter_cor_x1 == STN_267_x:
STN_v1 = 267
elif inter_cor_x1 == STN_277_x:
STN_v1 = 277
else:
print("no matching stations in v1")
if inter_cor_x2 == STN_229_x:
STN_v2 = 229
elif inter_cor_x2 == STN_235_x:
STN_v2 = 235
elif inter_cor_x2 == STN_242_x:
STN_v2 = 242
elif inter_cor_x2 == STN_251_x:
STN_v2 = 251
elif inter_cor_x2 == STN_267_x:
STN_v2 = 267
elif inter_cor_x2 == STN_277_x:
STN_v2 = 277
else:
print("no matching stations in v2")
if inter_cor_x3 == STN_229_x:
STN_v3 = 229
elif inter_cor_x3 == STN_235_x:
STN_v3 = 235
elif inter_cor_x3 == STN_242_x:
STN_v3 = 242
elif inter_cor_x3 == STN_251_x:
STN_v3 = 251
elif inter_cor_x3 == STN_267_x:
STN_v3 = 267
elif inter_cor_x3 == STN_277_x:
STN_v3 = 277
else:
print("no matching stations in v1")
inter_wind_v1 = wind[wind['STN'] == STN_v1]
inter_wind_v2 = wind[wind['STN'] == STN_v2]
inter_wind_v3 = wind[wind['STN'] == STN_v3]
df_add = inter_wind_v1.FF.add(inter_wind_v1.FF, fill_value=0)
inter_wind_v1 = wind[wind['STN'] == STN_v1].reset_index()
# =============================================================================
# print(f"inter_wind_v1:\n{inter_wind_v1}")
# =============================================================================
inter_wind_v2 = wind[wind['STN'] == STN_v2].reset_index()
# =============================================================================
# print(f"inter_wind_v2:\n{inter_wind_v2}")
# =============================================================================
inter_wind_v3 = wind[wind['STN'] == STN_v3].reset_index()
# =============================================================================
# print(f"inter_wind_v3:\n{inter_wind_v3}")
# =============================================================================
interp = pd.DataFrame(inter_wind_v1['FF']).rename(columns = {'FF':'inter_wind_v1'})
interp['inter_wind_v2'] = inter_wind_v2['FF']
interp['inter_wind_v3'] = inter_wind_v3['FF']
# =============================================================================
# print(interp)
# =============================================================================
def doInterp(inter_wind_v1, inter_wind_v2, inter_wind_v3):
return bary_interpol (inter_cor_x1, inter_cor_x2, inter_cor_x3, inter_cor_y1,
inter_cor_y2, inter_cor_y3, inter_cor_px, inter_cor_py,
inter_wind_v1, inter_wind_v2, inter_wind_v3)
interp['inter_wind_pv'] = interp.apply(lambda x: doInterp(x['inter_wind_v1'], x['inter_wind_v2'], x['inter_wind_v3']), axis=1)
print(interp)
filepath = Path('C:/Users/###/sensoren/data/out.csv')
filepath.parent.mkdir(parents=True, exist_ok=True)
interp.to_csv(filepath)
os.makedirs('C:/Users/###/sensoren/data/', exist_ok=True)
interp.to_csv('C:/Users/###/sensoren/data/out2.csv')
data_s = np.array([
[RD_229],
[RD_235],
[RD_242],
[RD_251],
[RD_267],
[RD_277]
])
data_p = np.array([
[P]
])
data_i = np.array([
[ inter_cor_x1, inter_cor_y1],
[ inter_cor_x2, inter_cor_y2],
[ inter_cor_x3, inter_cor_y3]
])
x_s, y_s = data_s.T
x_p, y_p = data_p. T
x_i, y_i = data_i. T
plt.scatter( x_p, y_p, c="red")
plt.scatter( x_s, y_s, c="blue")
plt.scatter( x_i, y_i, c="green")
plt.show()
## COPY POINT 3 ##
P = transformer.transform(50.00, 5.00) #location outside any triangle
## find best triangle or point
## find best triangle or point
triangle_fit = []
for triangle in itertools.combinations(points, 3):
p1, p2, p3 = triangle
if is_inside(*p1, *p2, *p3, *P):
triangle_fit.append(triangle)
# best perimeter
triangle_fit.sort(key=perimeter)
triangle_fit.append(triangle)
triangle_fit.sort(key=perimeter)
inter_cor = triangle_fit[0]
inter_cor_x1 = inter_cor[0][0]
inter_cor_x2 = inter_cor[1][0]
inter_cor_x3 = inter_cor[2][0]
inter_cor_y1 = inter_cor[0][1]
inter_cor_y2 = inter_cor[1][1]
inter_cor_y3 = inter_cor[2][1]
inter_cor_px = P[0]
inter_cor_py = P[1]
print(inter_cor_x1)
print(inter_cor_x2)
print(inter_cor_x3)
print(inter_cor_y1)
print(inter_cor_y2)
print(inter_cor_y3)
if inter_cor_x1 == STN_229_x:
STN_v1 = 229
elif inter_cor_x1 == STN_235_x:
STN_v1 = 235
elif inter_cor_x1 == STN_242_x:
STN_v1 = 242
elif inter_cor_x1 == STN_251_x:
STN_v1 = 251
elif inter_cor_x1 == STN_267_x:
STN_v1 = 267
elif inter_cor_x1 == STN_277_x:
STN_v1 = 277
else:
print("no matching stations in v1")
if inter_cor_x2 == STN_229_x:
STN_v2 = 229
elif inter_cor_x2 == STN_235_x:
STN_v2 = 235
elif inter_cor_x2 == STN_242_x:
STN_v2 = 242
elif inter_cor_x2 == STN_251_x:
STN_v2 = 251
elif inter_cor_x2 == STN_267_x:
STN_v2 = 267
elif inter_cor_x2 == STN_277_x:
STN_v2 = 277
else:
print("no matching stations in v2")
if inter_cor_x3 == STN_229_x:
STN_v3 = 229
elif inter_cor_x3 == STN_235_x:
STN_v3 = 235
elif inter_cor_x3 == STN_242_x:
STN_v3 = 242
elif inter_cor_x3 == STN_251_x:
STN_v3 = 251
elif inter_cor_x3 == STN_267_x:
STN_v3 = 267
elif inter_cor_x3 == STN_277_x:
STN_v3 = 277
else:
print("no matching stations in v1")
inter_wind_v1 = wind[wind['STN'] == STN_v1]
inter_wind_v2 = wind[wind['STN'] == STN_v2]
inter_wind_v3 = wind[wind['STN'] == STN_v3]
df_add = inter_wind_v1.FF.add(inter_wind_v1.FF, fill_value=0)
inter_wind_v1 = wind[wind['STN'] == STN_v1].reset_index()
# =============================================================================
# print(f"inter_wind_v1:\n{inter_wind_v1}")
# =============================================================================
inter_wind_v2 = wind[wind['STN'] == STN_v2].reset_index()
# =============================================================================
# print(f"inter_wind_v2:\n{inter_wind_v2}")
# =============================================================================
inter_wind_v3 = wind[wind['STN'] == STN_v3].reset_index()
# =============================================================================
# print(f"inter_wind_v3:\n{inter_wind_v3}")
# =============================================================================
interp = pd.DataFrame(inter_wind_v1['FF']).rename(columns = {'FF':'inter_wind_v1'})
interp['inter_wind_v2'] = inter_wind_v2['FF']
interp['inter_wind_v3'] = inter_wind_v3['FF']
# =============================================================================
# print(interp)
# =============================================================================
def doInterp(inter_wind_v1, inter_wind_v2, inter_wind_v3):
return bary_interpol (inter_cor_x1, inter_cor_x2, inter_cor_x3, inter_cor_y1,
inter_cor_y2, inter_cor_y3, inter_cor_px, inter_cor_py,
inter_wind_v1, inter_wind_v2, inter_wind_v3)
interp['inter_wind_pv'] = interp.apply(lambda x: doInterp(x['inter_wind_v1'], x['inter_wind_v2'], x['inter_wind_v3']), axis=1)
print(interp)
filepath = Path('C:/Users/###/sensoren/data/out.csv')
filepath.parent.mkdir(parents=True, exist_ok=True)
interp.to_csv(filepath)
os.makedirs('C:/Users/###/sensoren/data/', exist_ok=True)
interp.to_csv('C:/Users/###/sensoren/data/out2.csv')
data_s = np.array([
[RD_229],
[RD_235],
[RD_242],
[RD_251],
[RD_267],
[RD_277]
])
data_p = np.array([
[P]
])
data_i = np.array([
[ inter_cor_x1, inter_cor_y1],
[ inter_cor_x2, inter_cor_y2],
[ inter_cor_x3, inter_cor_y3]
])
x_s, y_s = data_s.T
x_p, y_p = data_p. T
x_i, y_i = data_i. T
plt.scatter( x_p, y_p, c="red")
plt.scatter( x_s, y_s, c="blue")
plt.scatter( x_i, y_i, c="green")
plt.show()
There are 2 repetitions in this script to see if it works. It creates the following 3 plots:
plotif blockquote still dont work: https://imgur.com/3327YX3
And it gives the following dataframe output:
2.3863918158
inter_wind_v1 inter_wind_v2 inter_wind_v3 inter_wind_pv
0 20.0 80.0 30.0 42.620259
1 30.0 60.0 40.0 42.816514
2 20.0 60.0 40.0 39.097197
3 20.0 40.0 40.0 32.561367
4 20.0 50.0 30.0 32.816514
... ... ... ...
10459 40.0 80.0 50.0 56.084429
10460 50.0 80.0 50.0 59.803745
10461 50.0 90.0 60.0 66.084429
10462 50.0 80.0 60.0 62.816514
10463 60.0 80.0 60.0 66.535830
[10464 rows x 4 columns]
inter_wind_v1 inter_wind_v2 inter_wind_v3 inter_wind_pv
0 40.0 30.0 30.0 37.356121
1 30.0 40.0 20.0 31.560683
2 20.0 40.0 30.0 24.746160
3 30.0 40.0 30.0 32.102281
4 20.0 30.0 30.0 22.643879
... ... ... ...
10459 50.0 50.0 50.0 50.000000
10460 50.0 50.0 50.0 50.000000
10461 40.0 60.0 50.0 44.746160
10462 60.0 60.0 60.0 60.000000
10463 60.0 60.0 50.0 59.458402
[10464 rows x 4 columns]
152273.8331310129
154737.80483951585
208796.2620382847
600780.6766655241
545552.3863918158
602766.4003272419
inter_wind_v1 inter_wind_v2 inter_wind_v3 inter_wind_pv
0 40.0 30.0 30.0 -20.684217
1 30.0 40.0 20.0 105.480717
2 20.0 40.0 30.0 148.766684
3 30.0 40.0 30.0 98.082467
4 20.0 30.0 30.0 80.684217
... ... ... ...
10459 50.0 50.0 50.0 50.000000
10460 50.0 50.0 50.0 50.000000
10461 40.0 60.0 50.0 168.766684
10462 60.0 60.0 60.0 60.000000
10463 60.0 60.0 50.0 67.398250
[10464 rows x 4 columns]
As you can see the 2 scatter plot show the wanted results. The used stations change depending on (P). But when point P is outside of the triangles I would like to not interpolate but take the wind(FF) value of the station closest to point (P) It could be a code that look like this:
inter_cor=(min_distance(P[0], P[1], points))
But I am not able to implement this. I wonder if somewan knows how I can implement this. And that the dataframe output will be from just one station instead of 3.
I hope somewan can help me with this! I am almost halfway
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|