[Python] scipy curve_fit fails on specified function

I am trying to find the optimal projection parameters when transforming spherical coordinates into Cartesian ones for a given projection. There are four unknown projection parameters, in addition, the coordinates on the plane can change using the Helmert transformation, so I need to find 8 parameters for the function I have specified. An array of spherical coordinates and the corresponding coordinates on the plane are known (12 points). I decided to use curve_fit from the scipy package for this purpose. However, this function displays inaccurate results. Below is the script code and the program output. Optimal parameter values in code comments

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)


Output:

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


as you can see the results are inaccurate because pcov has big values and angles have 7 degrees deviation. I will be grateful for any help, also if there are recommendations to use other means, can you tell about them.

I tried scaling variables, changing the power of 10, changing sigma, decreasing/increasing the number of points, removing the YTranslate parameter

Continue reading...