Page 1 of 1
Math related question, level 2
Posted: Tue Dec 08, 2009 11:33 pm
by 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
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 12:28 am
by Demivec
Michael Vogel wrote: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) .
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 2:54 am
by 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)
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 7:32 am
by Michael Vogel
Demivec wrote: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...
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 8:10 am
by 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
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 11:59 am
by Michael Vogel
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
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 1:10 pm
by 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!
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 4:18 pm
by Michael Vogel
dioxin wrote: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
Re: Math related question, level 2
Posted: Wed Dec 09, 2009 8:31 pm
by 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
Re: Math related question, level 2
Posted: Fri Dec 11, 2009 3:58 pm
by 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)
Re: Math related question, level 2
Posted: Tue Dec 15, 2009 9:33 am
by 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