Scipy curve_fit schlägt bei der angegebenen Funktion fehl ⇐ Python

[Anonymous]

Ich versuche, die optimalen Projektionsparameter zu finden, wenn die kugelförmigen Koordinaten für eine bestimmte Projektion in kartesische Umstände umgewandelt werden. Die mathematische Beschreibung der Projektion wird in den Bildern dargestellt (Radius wird aufgrund der Helmert -Transformation weggelassen). Parameter bereiten die Koordinatenberechnung vor. Es gibt vier unbekannte Projektionsparameter. Darüber hinaus können sich die Koordinaten in der Ebene unter Verwendung der Helmert -Transformation ändern. Daher muss ich 8 Parameter für die von mir angegebene Funktion finden. Eine Reihe kugelförmiger Koordinaten und die entsprechenden Koordinaten in der Ebene sind (12 Punkte) bekannt. Ich habe mich für diesen Zweck entschieden, Curve_Fit aus dem Scipy -Paket zu verwenden. Diese Funktion zeigt jedoch ungenaue Ergebnisse an. Unten finden Sie den Skriptcode und die Programmausgabe. Optimale Parameterwerte in Code-Kommentaren < /p>
import numpy as np
from scipy.optimize import curve_fit

def geographic2ProjAndHelmert(VRadianCoord, lat_0, lon_0, lat_st_1, lat_st_2, a, b, TranslationX, TranslationY):
#VRadianCoord - [lon1 (x), lat1 (y), lon2, lat2, ...]; a - scale * cos(φ); b - scale * sin(φ)

VRadianCoord=VRadianCoord.reshape((-1, 2))

#scale variables
a = a * 10
b = b * 10
lat_0 = lat_0 / 1000
lon_0 = lon_0 / 1000
lat_st_1 = lat_st_1 / 1000
lat_st_2 = lat_st_2 / 1000

#Calculate Cartesian
cos1 = np.cos(lat_st_1)
cos2 = np.cos(lat_st_2)
CompositeTanCalc0 = np.tan(np.pi/4.0 + lat_0/2.0)
CompositeTanCalc1 = np.tan(np.pi/4.0 + lat_st_1/2.0)
CompositeTanCalc2 = np.tan(np.pi/4.0 + lat_st_2/2.0)
n = (np.log(cos1/cos2)) / (np.log(CompositeTanCalc2/CompositeTanCalc1))
F = cos1*CompositeTanCalc1**n/n
r = F/np.tan(np.pi/4.0 + VRadianCoord[:,1]/2.0)**n
r_0 = F/CompositeTanCalc0**n
ProjCoord = np.array([r*np.sin(n*(VRadianCoord[:,0]-lon_0)), r_0-r*np.cos(n*(VRadianCoord[:,0]-lon_0))])

#Calculate Helmert Transformation
RotationAndScale = np.array([[a, -b], [b, a]])
EndCoord = np.copy(np.dot(RotationAndScale, ProjCoord).T)
EndCoord[:,0] = np.add(EndCoord[:,0], TranslationX)
EndCoord[:,1] = np.add(EndCoord[:,1], TranslationY)

return EndCoord.ravel()

RadCoord = np.array([0.73303829, 1.15191731, 1.01229097, 1.13446401, 0.9250245, 0.97738438, 0.71558499, 1.01229097, 0.87266463, 1.15191731, 0.83775804, 0.99483767, 0.95993109, 1.06465084, 0.71558499, 1.08210414, 0.87266463, 1.06465084, 0.78539816, 1.02974426, 0.80285146, 1.08210414, 0.89011792, 1.01229097])
#format [x1, y1, x2, y2, ...]

Ydata = np.array([897.66512043, -4893.42010515, 585.79116218, -470.03157748, 6683.05857945, -272.07474213, 6149.61630033, -4694.83472251, 550.95105566, -2763.00971689, 6453.46869002, -2176.47012956, 3359.83097409, -508.82821497, 3546.62200096, -4926.08133397, 3761.31345969, -2047.13507678, 5341.79985183, -3415.91999683, 3355.96753336, -3399.9554403, 5595.98729439, -1281.99094122])

p0 = np.radians(55)*1000, np.radians(85)*1000, np.radians(40)*1000, np.radians(69)*1000, -3705.8, 4627, -2887, 7500
#p0_optimal = np.radians(60)*1000, np.radians(90)*1000, np.radians(48)*1000, np.radians(64)*1000, -2705.8, 2627, -1887, 9500

popt, pcov = curve_fit(geographic2ProjAndHelmert, RadCoord, Ydata, p0, maxfev=100000, bounds=([np.radians(30)*1000,np.radians(60)*1000,np.radians(30)*1000, np.radians(40)*1000, -6000, -6000, -6000, -10000],[np.radians(70)*1000,np.radians(100)*1000,np.radians(70)*1000, np.radians(70)*1000, 6000, 6000, 6000, 10000]))
popt[0:4] = np.degrees(popt[0:4])/1000
print("popt = ", popt, "\n")
print("pcov = ", pcov)
< /code>
Ausgabe: < /p>
popt = [ 59.20915123 83.65741635 40.57766132 69.72091072
-2508.47749214 2936.1387793 -73.51463145 8172.85715372]

pcov = [[ 9.18230672e+12 -2.08906927e+13 1.99021988e+12 -1.51118796e+12
5.20958799e+13 4.21175110e+13 -1.30500456e+12 -5.27326558e+14]
[-2.08906927e+13 4.76711630e+13 -5.13324974e+12 3.89771178e+12
-1.19230289e+14 -9.56985287e+13 1.19399847e+12 1.20179982e+15]
[ 1.99021988e+12 -5.13324974e+12 6.92449865e+12 -5.25781389e+12
1.66182697e+13 5.88100103e+12 7.24698818e+12 -1.23109006e+14]
[-1.51118796e+12 3.89771178e+12 -5.25781389e+12 3.99229002e+12
-1.26183630e+13 -4.46548828e+12 -5.50267812e+12 9.34775221e+13]
[ 5.20958799e+13 -1.19230289e+14 1.66182697e+13 -1.26183630e+13
3.00448606e+14 2.36727395e+14 1.38861694e+12 -3.00208240e+15]
[ 4.21175110e+13 -9.56985287e+13 5.88100103e+12 -4.46548828e+12
2.36727395e+14 1.95183572e+14 -7.51770717e+12 -2.41695451e+15]
[-1.30500456e+12 1.19399847e+12 7.24698818e+12 -5.50267812e+12
1.38861694e+12 -7.51770717e+12 2.22666466e+13 4.90988220e+13]
[-5.27326558e+14 1.20179982e+15 -1.23109006e+14 9.34775221e+13
-3.00208240e+15 -2.41695451e+15 4.90988220e+13 3.03138515e+16]]
< /code>
Wie Sie sehen können, sind die Ergebnisse ungenau, da PCOV große Werte und Winkel auf 7 Grad -Abweichungen aufweist. Ich werde für jede Hilfe dankbar sein, auch wenn es Empfehlungen gibt, die andere Mittel verwenden, können Sie über sie erzählen. parameters by a special preliminary coordinate transformation:

[*]translation to the center of mass

[*]scaling (the largest distance from one of the points to the center of coordinates is 1)

rotation (the point that is Am weitesten aus dem Zentrum der Koordinaten liegt auf einer der Achsen).>