経緯線交点のピクセル座標を読み取った結果 （左から、経度、緯度、X座標、下向きY座標）
$ cat aupq78.txt
#proj +proj=stere +lat_ts=60 +lat_0=90 +lon_0=50 +k_0=1.0 +x_0=-1962800 +y_0=7405800 +ellps=GRS80
80/30/384/391
80/40/685/217
100/30/789/872
100/10/261/1501
110/60/1531/213
110/10/636/1765
110/50/1381/475
120/10/1052/1959
130/60/1806/313
130/50/1753/611
130/40/1697/926
140/60/1952/327
140/40/1952/949
140/30/1951/1295
140/20/1952/1682
160/60/2239/277
160/40/2452/860
160/30/2571/1186
160/20/2702/1549
170/60/2371/213
170/50/2522/475
170/40/2682/752
180/50/2685/360
ここから当てはめを行うスクリプト
$ cat fitmap.rb
proj = nil
sigma = 2.0
pts = []
for line in ARGF
case line
when /^#proj\s+/ then
proj = $'
next
when /^#sigma\s+/ then
sigma = $'.to_f
next
when /^#/ then next
end
lon, lat, x, y = line.strip.split(/[\/\s,]+/).map{|s| s.to_f}
pxy = (`echo #{lon} #{lat} | proj #{proj}`)
px, py = pxy.strip.split(/\s+/).map{|s| s.to_f}
pt = {:lon => lon, :lat => lat, :x => x, :y => y, :px => px, :py => py}
pts.push pt
end
def stat(pts, zx, zy, sig)
ss = sx = sy = 0.0
for pt in pts
wt = (sig ** -2)
ss += wt
sx += pt[zx] * wt
sy += pt[zy] * wt
end
sxoss = sx / ss
st2 = b = 0.0
for pt in pts
t = (pt[zx] - sxoss) / sig
st2 += t * t
b += t * pt[zy] / sig
end
b /= st2
a = (sy - sx * b) / ss
siga = Math.sqrt((1.0 + sx * sx / (ss * st2)) / ss)
sigb = Math.sqrt(1.0 / st2)
chi2 = 0.0
for pt in pts
e2 = ((pt[zy] - a - b * pt[zx]) / sig) ** 2
STDERR.puts [e2, pt].inspect if e2 > 4.0
chi2 += e2
end
{:a => a, :b => b, :siga => siga, :sigb => sigb, :chi2 => chi2}
end
def repstat h
printf("a:%-10.1f (%-10.2f) b:%13.4e (%13.3e) a/b:%-8.0f 1/b:%-6.0f chi2:%g\n",
h[:a], h[:siga], h[:b], h[:sigb], h[:a]/h[:b], 1.0/h[:b], h[:chi2])
end
repstat stat(pts, :px, :y, sigma)
repstat stat(pts, :py, :x, sigma)
内容は煎じ詰めて言うと、縦横方向バラバラに線形回帰をしていて、各点の測定精度は標準偏差2ピクセルとの仮定をしています。 極のピクセル番号 a とスケール（1ピクセルの地上長の逆数）bについてそれぞれ推定の標準偏差が示されます。 当てはめ精度が悪くて 2 sigma 以上のところは点の情報が表示されます。 stat 関数は Numerical Recipes in C のアルゴリズムでピクセル測定精度を一定と仮定して Ruby で書き直したものです
実行してみるとこうなるわけです
$ ruby fitmap.rb aupq78.txt
a:-0.0 (0.76 ) b: 2.6438e-04 ( 2.014e-07) a/b:-27 1/b:3782 chi2:10.1364
a:-0.0 (1.04 ) b: 2.6347e-04 ( 1.447e-07) a/b:-51 1/b:3795 chi2:7.72033
つまり X 軸のピクセル間隔は 3782 m, Y 軸は 3796 m になります。
もういっこいきましょう。
$ cat auas50.txt
#proj +proj=stere +lat_ts=60 +lat_0=90 +lon_0=50 +k_0=1.0 +x_0=-1696000 +y_0=7268100 +ellps=GRS80
70/40/389/41
80/30/251/332
80/40/466/206
80/50/666/90
90/20/167/696
90/30/379/517
90/40/571/357
90/50/748/208
90/60/914/69
100/10/161/1129
100/20/362/890
100/30/540/677
100/40/702/486
100/50/849/309
100/60/988/143
110/0/249/1633
110/10/431/1319
110/20/588/1048
110/30/726/807
110/40/851/591
110/50/966/391
110/60/1074/204
110/70/1177/26
120/10/729/1457
120/20/836/1163
120/30/931/903
120/40/1016/668
120/50/1096/452
120/60/1170/248
120/70/1240/53
130/10/1048/1543
130/20/1101/1235
130/30/1150/962
130/40/1194/714
130/50/1234/489
130/60/1271/275
130/70/1307/72
140/10/1376/1572
140/20/1377/1259
140/30/1376/981
140/40/1377/729
140/50/1376/502
140/60/1376/285
140/70/1376/79
150/10/1704/1542
150/20/1650/1234
150/30/1602/961
150/40/1558/715
160/10/2023/1458
160/20/1916/1163
160/30/1821/903
160/50/1657/451
160/60/1583/249
160/70/1511/54
170/30/2027/808
170/40/1901/591
170/50/1786/391
170/60/1678/204
170/70/1575/25
180/40/2051/485
180/50/1904/309
180/60/1764/143
190/50/2004/209
190/60/1838/70
200/50/2087/90
実行結果はこうなります
$ ruby fitmap.rb auas50.txt
a:-0.0 (0.41 ) b: 1.8992e-04 ( 1.005e-07) a/b:-113 1/b:5265 chi2:27.1943
[4.19913391375125, {:x=>249.0, :y=>1633.0, :lon=>110.0, :px=>8577371.02, :lat=>0.0, :py=>1336766.48}]
a:0.0 (0.60 ) b: 1.8932e-04 ( 8.806e-08) a/b:124 1/b:5282 chi2:19.095
ひとつエラーになってますが、これだけ点数があると、パラメタ推定の標準偏差がえらく小さくなってしまうので、そこからすこしでもずれると大きく評価されていまうようです。