Math related question, level 2

[Michael Vogel]

Given is a rectangle map on which two points are marked and the coordinates of these two points are known, how to calculate the coordinates of the corners and the orientation of the map (direction of to the north pole)?

+---------------------------------+
|·································|
|·······A························| A: N 48°30' | W 5°20'
|·································|
|·································|
|·································|
|··················B·············| B: N 48°00' | W 6°20'
|·································|
|·································|
+---------------------------------+

Seems not so easy

[Demivec]

Given is a rectangle map on which two points are marked and the coordinates of these two points are known, how to calculate the coordinates of the corners and the orientation of the map (direction of to the north pole)?

Not enough information is given to solve the problem. What is needed: size of map and coordinate of one corner, and 2 more points that are parallel to one edge of the map (one of the points may be the already included corner). None of the additional points would be needed if the coordinates of points A and B was given also in absolute map coordinates(x,y) in addition to (latitude,longitude) .

[citystate]

the orientation of the map is probably the only thing you could answer (as A is to the NE of B according to the coordinates given, North should be towards the left of the map)

[Michael Vogel]

Not enough information is given to solve the problem. What is needed: size of map and coordinate of one corner, and 2 more points that are parallel to one edge of the map (one of the points may be the already included corner). None of the additional points would be needed if the coordinates of points A and B was given also in absolute map coordinates(x,y) in addition to (latitude,longitude) .

Sorry, everytime I start translating the text from my mother language into english, a lot of information get lost

The map itself is a bitmap (e.g. a satellite photo or a scan of a real map), so you have a pixel width and height of the whole "map" and the two points are also given by "relative pixel coordinates" within that bitmap image.

The "world coordinates" are just from our earth.

Michael

PS maybe it is clearer when I say I want to create "garmin custom maps" which allow to show bitmaps in a GPS device - I started yesterday with the main program (see here) but failed to do a correct calculation of the coordinates...

[idle]

you need to know what the projection type is, is it a Mercator projection or some other beast?\

google maps use
http://en.wikipedia.org/wiki/Mercator_projection

google satelite use
http://en.wikipedia.org/wiki/Plate_carr ... projection

[Michael Vogel]

you need to know what the projection type is, is it a Mercator projection or some other beast?\

google maps use
http://en.wikipedia.org/wiki/Mercator_projection

google satelite use
http://en.wikipedia.org/wiki/Plate_carr ... projection

You're right, if the map is large enough, also these things will influence the whole thing (I had to do such things already for my Forerunner tool to get google maps displayed for specific areas)...

...but for smaller maps (if not to near to the poles) with a size of few kilometers this (hopefully) could be ignored, so I want to start with that using even non spherical geometry...

I thought that the first step would be to find the point of intersection of two circles. The first one has its center in point A an the radius is the delta of the latitudes and the second has its center in B and the radius is the delta of the longitudes. But this won't help me to get the orientation of the map. As soon this is found, the rest (at least in 2D) should be easy.

If these results are not precise enough, everything will go even worse

Michael

[dioxin]

A and B are close enough together that we can ignore the curvature of the Earth and still get a good result so just use basic trigonometry.

A is North of B by 48d30m - 48d00m = 0d30m = 0.5d
A is East of B by 6d20m - 5d20m = 1d00m = 1.0d

North is therefore off to the bottom left of the grid at an angle of atan(1.0/0.5) = atan(2) = 63.435d to the line AB.


Assuming the dots/lines on the grid represent equal distances within the grid.
Using Pythagoras:
In grid units AB is length sqr(4^2 + 11^2) = 11.704
In degrees AB is length sqr(1^2 + 0.5^2) = 1.118d
Each degree is then 11.704 / 1.118 = 10.469 grid units.

Plot a point D vertically down from A and horizontally left from B.
The angle ABD is then atan(4/11) = 19.983d

North is then 63.435d - 19.983d = 43.452d to the grid, i.e angle DBN = 43.452d

Now, in turn, take each corner of the grid and using it, form a right angled triangle with A and the line heading North from A.
That triangle has a known hypotenuse (Pythagoras) and a known angle (atan on the angle plus the angle DBN) so the distance of the corner North and West from A can easily be found.

I tried posting this yesterday but the forum kept saying I have no permission to access this forum!

[Michael Vogel]

A and B [...]
I tried posting this yesterday but the forum kept saying I have no permission to access this forum!

Thanks, but now you have to wait another day - it will take a while to check (and understand :roll:) that

Michael

[Michael Vogel]

Still have problems to get the right results, maybe there is an issue because of the circles of latitudes get shorter north and sorth from the equator?

Code: Select all

For grad=0 To 90 Step 10
	Debug Cos(grad/180*#PI)*40000
Next grad

New Target is to get the same results like a program I've just found in the internet, which calculated the following (hopefully correct) results based on a bitmap with the height and wisth of 512 (from 0|0 to 512|512):

Point A is on 99|99 and was defined as 48°N 20°E
Point B is on 349|299 and was defined as 46°N 25°E

And here are the results of the program:
0|0: 49.5279075155°N 18.5911840832°E
511|511: 43.1910710990°N 27.8537203726°E
Rotation: 8.3458831017

This values have been seen also on the output, but are not important for me:
Aspect ratio: 1.46169723829597
D: 1375.37040680529 metres/pixel
T: 0.145663139110928

[Michael Vogel]

Dioxin brought me near to the solution, but I still fail to get the correct results

Here's what I have so far...

Code: Select all

Macro Grad(x)
	(x)*180/#PI
EndMacro
Macro DebugN
	Debug ""
	Debug "## "+
EndMacro

; 0|0: 49.5279075155°N 18.5911840832°E
; 511|511: 43.1910710990°N 27.8537203726°E
; Rotation: 8.3458831017
; Aspect ratio: 1.46169723829597
; D: 1375.37040680529 metres/pixel
; T: 0.145663139110928

Global Dim Lat.d(2)
Global Dim Lon.d(2)
Global Dim X(2)
Global Dim Y(2)

Read.l X(0)
Read.l Y(0)

For i=1 To 2
	Read.l X(i)
	Read.l Y(i)
	Read.d Lat(i)
	Read.d Lon(i)
Next i

DataSection
	Data.l 512,512
	Data.l 99,99
	Data.d 48.0
	Data.d 20.0
	Data.l 349,299
	Data.d 46.0
	Data.d 25.0
EndDataSection

DataSection
	Data.l 512,512
	Data.l 23,30
	Data.d 3.0
	Data.d 7.0
	Data.l 67,60
	Data.d 10.0
	Data.d 15.0
EndDataSection


; Deltas
Lat(0)=Lat(2)-Lat(1)
Lon(0)=Lon(2)-Lon(1)

X(0)=X(2)-X(1)
Y(0)=Y(2)-Y(1)

DebugN "Delta (World):"
Debug lat(0)
Debug lon(0)

; Distances:
Wd.d=Sqr(Lat(0)*Lat(0)+Lon(0)*Lon(0))
Md.d=Sqr(X(0)*X(0)+Y(0)*Y(0))

DebugN "Distance: "
Debug Wd
Debug Md

; Angles
Ma.d=ASin(X(0)/Md)
Mb.d=ASin(Y(0)/Md)
Wa.d=ASin(Lat(0)/Wd)
Wb.d=ASin(Lon(0)/Wd)
Dif1.d=Wa-Ma
Dif2.d=Wb-Mb

DebugN "Angles: "
Debug Grad(Ma)
Debug Grad(Mb)
Debug Grad(Wa); Nord/Süd
Debug Grad(Wb); West/Ost

DebugN "Check Coordinates:"
Debug Cos(Wb)*wd
Debug Cos(Wa)*wd

; Calculate 0|0:
x.d=-X(1)
y.d=-Y(1)
r.d=Sqr(x*x+y*y)

DebugN "Distance to 0|0:"
Debug r

; Direction to 0|0:
NeuA=ASin(x/r)
NeuA-Wa

DebugN "0|0 World Coordinates:"
Debug Sin(NeuA)*r*wd/md+Lat(1)
Debug Cos(NeuA)*r*wd/md+Lon(1)

[Michael Vogel]

Got one step closer again - the main problem is to get the correct rotation of the map, and for this the ratio of the x-scale to y-scale is needed...

But this should be possible now with this procedure which calculates the distance between two points on the earth surface:

Code: Select all

Procedure.d Distance(Lat1.d,Lon1.d,Lat2.d,Lon2.d)

	Protected a.d=6378137;							Große Achse des Ellipsoids
	Protected b.d=6356752.3142;					Kleine Achse des Ellipsoids
	Protected f.d=1/298.257223563;					Abflachung des WGS-84 Ellipsiods=(A-B)/A

	Protected L.d=Radiant(Lon2-Lon1);				Differenz der Längengrade

	Protected U1.d=ATan((1-f)*Tan(Radiant(Lat1)));	Geodätische Länge
	Protected U2.d=ATan((1-f)*Tan(Radiant(Lat2)));	Geodätische Länge

	Protected SinU1.d=Sin(U1)
	Protected CosU1.d=Cos(U1)
	Protected SinU2.d=Sin(U2)
	Protected CosU2.d=Cos(U2)

	Protected Lambda.d=L;							erste Näherung
	Protected LambdaP.d;							letzter Näherungswert
	Protected IterLimit=100;							Anzahl der maximalen "Annäherungsversuche"

	#Epsilon=1E-12;									gewünschte Genauigkeit

	Repeat

		Protected SinLambda.d=Sin(Lambda)
		Protected CosLambda.d=Cos(Lambda)

		Protected U2Lambda.d=CosU2*SinLambda
		Protected U1U2Lambda.d=CosU1*SinU2-SinU1*CosU2*CosLambda

		Protected SinSigma.d=U2Lambda*U2Lambda+U1U2Lambda*U1U2Lambda
		If SinSigma<=0
			ProcedureReturn -1;						Co-incident Points
		Else
			SinSigma=Sqr(SinSigma)
		EndIf
		;Debug "ok SinSigma "+StrD(SinSigma,20)


		Protected CosSigma.d=SinU1*SinU2+CosU1*CosU2*CosLambda;
		Protected Sigma.d=ATan2(SinSigma,CosSigma)


		Protected SinAlpha.d=CosU1*CosU2*SinLambda/SinSigma;

		Protected CosSqAlpha.d=1-SinAlpha*SinAlpha;

		Protected Cos2SigmaM.d=0; 					Equatorial line
		If CosSqAlpha<>0
			Cos2SigmaM.d=CosSigma-2*SinU1*SinU2/CosSqAlpha;
		EndIf

		Protected C.d=f/16*CosSqAlpha*(4+f*(4-3*CosSqAlpha));

		LambdaP=Lambda;
		Lambda=L+(1-C)*f*SinAlpha*(sigma + C*SinSigma*(Cos2SigmaM+C*CosSigma*(-1+2*Cos2SigmaM*Cos2SigmaM)))

	Until (IterLimit=0) Or (Abs(Lambda-LambdaP)<#Epsilon)

	If iterLimit=0
		ProcedureReturn 0;NaN  // formula failed To converge
	EndIf

	Protected uSq.d=CosSqAlpha*(a*a-b*b)/(b*b);

	Protected A_.d=1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
	Protected B_.d=uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
	Protected deltaSigma.d=B_*SinSigma*(Cos2SigmaM+B_/4*(CosSigma*(-1+2*Cos2SigmaM*Cos2SigmaM)- B_/6*Cos2SigmaM*(-3+4*SinSigma*SinSigma)*(-3+4*Cos2SigmaM*Cos2SigmaM)));

	Protected s.d=b*A_*(sigma-deltaSigma);

	ProcedureReturn s

EndProcedure