Ahoj!
> > mam ;-). Bohuzel je v pascalu.
>
> Prepsat to do C++ nebo Javy by nemel byt problem :) Posli do mailing listu ...
Javu uz pry nekdo ma... (Nemelo by tohle umet proj?). Pascali zdrojaky
prilozeny.
> > V jakym formatu ma byt vystup? Predpokladam ze .osm soubor pro
> > josm... kde je we wikine popsany format adresniho bodu / muze mi nekdo
> > par adresnich bodu poslat? (Jsem na *hodne* pomaly lince :-().
>
> na git.wz.cz doslo misto, takze to dam jinam:
>
> http://nasa.kx.cz/doc.7z - dokumentace, 70kb
> http://nasa.kx.cz/demo_sql.7z - demo data v sql formatu, 2.5 mb
> http://nasa.kx.cz/demo.7z - demo data v csv formau, 2.5 mb
...dik, ty uz ma. Ja myslel adresni bod ze streetmapy abych vedel jak
to ma vypadat. Vim ze nekde v praze uz jsou, ale uz si nepamatuju kde
a nemuzu to najit :-(.
--
(english) http://www.livejournal.com/~pavelmachek
(cesky, pictures)
http://atrey.karlin.mff.cuni.cz/~pavel/picture/horses/blog.html
{
Copyright 2005 Zdenek Hrdina, distribute under GPLv2
}
procedure jtsk_wgs( X,Y,Hel:double; var B,L,H:double);
{Vypocet zemepisnych souradnic v systemu WGS-84 z rovinnych souradnic
S-JTSK a elipsoidicke vysky}
procedure transformace_BLH(var B,L,H: double);
{Transformace zemepisnych souradnic z JTSK do WGS}
var lat,lon,alt,x1,y1,z1,x2,y2,z2:double;
procedure transformace(xs,ys,zs:double; var xn,yn,zn:double);
{transformace pravouhlych souradnic}
const
{koeficienty transformace ze systemu S-JTSK do systemu WGS-84}
dx=570.69; dy=85.69; dz=462.84; {posunuti}
wz=-5.2611/3600*pi/180; wy=-1.58676/3600*pi/180; wx=-4.99821/3600*pi/180; {rotace}
m=3.543e-6; {meritko}
begin
xn:=dx+(1+m)*(xs+wz*ys-wy*zs);
yn:=dy+(1+m)*(-wz*xs+ys+wx*zs);
zn:=dz+(1+m)*(wy*xs-wx*ys+zs);
end;
procedure BLH_xyz(B,L,H:double; var x,y,z:double);
{vypocet pravouhlych souradnic ze zemepisnych souradnic}
const
{parametry Besselova elipsoidu}
a=6377397.15508; f_1=299.152812853;
var
ro,e2:double;
begin
e2:=1-sqr(1-1/f_1); ro:=a/sqrt(1-e2*sqr(sin(B)));
x:=(ro+H)*cos(B)*cos(L); y:=(ro+H)*cos(B)*sin(L); z:=((1-e2)*ro+H)*sin(B);
end;
procedure xyz_BLH(x,y,z:double; var B,L,H:double);
{vypocet zemepisnych souradnic z pravouhlych souradnic}
const
{parametry elipsoidu WGS-84}
a=6378137.0; f_1=298.257223563;
var
a_b,e2,theta,st,ct,p,t:double;
begin
a_b:=f_1/(f_1-1); p:=sqrt(sqr(x)+sqr(y)); e2:=1-sqr(1-1/f_1);
theta:=arctan(z*a_b/p);st:=sin(theta);ct:=cos(theta);
t:=(z+e2*a_b*a*sqr(st)*st)/(p-e2*a*sqr(ct)*ct);
B:=arctan(t);
H:=sqrt(1+sqr(t))*(p-a/sqrt(1+(1-e2)*sqr(t)));
L:=2*arctan(y/(p+x));
end;
begin
BLH_xyz(B,L,H,x1,y1,z1);
transformace(x1,y1,z1,x2,y2,z2);
xyz_BLH(x2,y2,z2,B,L,H);
end;
procedure XY_BL(X,Y:double; var B,L: double);
{Vypocet zemepisnych souradnic z rovinnych souradnic}
const
a=6377397.15508; e=0.081696831215303;
n=0.97992470462083; konst_u_ro=12310230.12797036;
sinUQ=0.863499969506341; cosUQ=0.504348889819882;
sinVQ=0.420215144586493; cosVQ=0.907424504992097;
alfa=1.000597498371542; k=1.003419163966575;
var
ro,epsilon,D,sinS,cosS,sinU,cosU,cosDV,sinDV,S,sinV,cosV,sinB,t,pom:double;
begin
ro:=sqrt(sqr(x)+sqr(y));
epsilon:=2*arctan(y/(ro+x));
D:=epsilon/n; S:=2*arctan(exp(1/n*ln(konst_u_ro/ro)))-pi/2;
sinS:=sin(S);cosS:=cos(S);
sinU:=sinUQ*sinS-cosUQ*cosS*cos(D);cosU:=sqrt(1-sqr(sinU));
sinDV:=sin(D)*cosS/cosU;cosDV:=sqrt(1-sqr(sinDV));
sinV:=sinVQ*cosDV-cosVQ*sinDV;cosV:=cosVQ*cosDV+sinVQ*sinDV;
L:=2*arctan(sinV/(1+cosV))/alfa;
t:=exp(2/alfa*ln((1+sinU)/cosU/k));
pom:=(t-1)/(t+1);
repeat
sinB:=pom;
pom:=t*exp(e*ln((1+e*sinB)/(1-e*sinB)));
pom:=(pom-1)/(pom+1);
until abs(pom-sinB)<1e-15;
B:=arctan(pom/sqrt(1-sqr(pom)));
end;
begin
XY_BL(X,Y,B,L); H:=Hel; transformace_BLH(B,L,H);
end; { jtsk_wgs }
var X,Y,Hel:double; var B,L,H:double;
begin
readln(x, y, Hel);
jtsk_wgs(x, y, hel, b, l, h);
writeln(b/(2*pi)*360 :10:7, ' ', l/(2*pi)*360 :10:7, h :10:3);
end.{
Copyright 2005 Zdenek Hrdina, distribute under GPLv2
Usage: ( echo lat; echo lon; echo alt ) | wgs2jtsk
lat and lon are in degrees.
}
program prepocet;
procedure transformace(xs,ys,zs:double; var xn,yn,zn:double);
{transformace pravouhlych souradnic}
const
{koeficienty transformace ze systemu WGS-84 do systemu S-JTSK}
dx=-574.4; dy=-119.4; dz=-421.6;
wz=2.5e-5;wy=3.825e-6;wx=3.162e-5;
m=-7.39e-6;
begin
xn:=dx+(1+m)*(xs+wz*ys-wy*zs);
yn:=dy+(1+m)*(-wz*xs+ys+wx*zs);
zn:=dz+(1+m)*(wy*xs-wx*ys+zs);
end;
procedure BLH_xyz(B,L,H:double; var x,y,z:double);
{vypocet pravouhlych souradnic z geodetickych souradnic}
const
{parametry elipsoidu WGS-84}
a=6378137.0;f_1=298.257223563;
var
ro,e2:double;
begin
e2:=1-sqr(1-1/f_1); ro:=a/sqrt(1-e2*sqr(sin(B)));
x:=(ro+H)*cos(B)*cos(L);
y:=(ro+H)*cos(B)*sin(L);
z:=((1-e2)*ro+H)*sin(B);
end;
procedure xyz_BLH(x,y,z:double; var B,L,H:double);
{vypocet geodetickych souradnic z pravouhlych souradnic}
const
{parametry Besselova elipsoidu}
a=6377397.15508; f_1=299.152812853;
var
a_b,e2,theta,st,ct,p,t:double;
begin
a_b:=f_1/(f_1-1); p:=sqrt(sqr(x)+sqr(y)); e2:=1-sqr(1-1/f_1);
theta:=arctan(z*a_b/p);st:=sin(theta);ct:=cos(theta);
t:=(z+e2*a_b*a*sqr(st)*st)/(p-e2*a*sqr(ct)*ct);
B:=arctan(t);
H:=sqrt(1+sqr(t))*(p-a/sqrt(1+(1-e2)*sqr(t)));
L:=2*arctan(y/(p+x));
end; { xyz_BLH }
function WGS84toBessel(var Latitude, Longitude, Altitude: Double): Boolean;
//------------------------------------------------------------------------------
var B, L, H, xl, yl, zl, x2, y2, z2 : double;
begin
B := Latitude;
L := Longitude;
H := Altitude;
BLH_xyz(B,L,H,xl,yl,zl);
transformace(xl,yl,zl,x2,y2,z2);
xyz_BLH(x2,y2,z2,B,L,H);
Latitude := B;
Longitude := L;
Altitude := H;
end;
function BesseltoJTSK(Latitude, Longitude: Double; var C_X, C_Y : double): Boolean;
// vraci JTSK souřadnice pro zadanou Severní sířku a východní délku
// není třeba odečítat Ferro
//******************************************************************************
const a = 6377397.15508;
e = 0.081696831215303;
n = 0.97992470462083;
rho_0 = 12310230.12797036;
sinUQ = 0.863499969506341;
cosUQ = 0.504348889819882;
sinVQ = 0.420215144586493;
cosVQ = 0.907424504992097;
alfa = 1.000597498371542;
k_2 = 1.00685001861538;
var rho, eps, V, t, B, L : Double;
sinS, cosS, sinU, cosU, cosDV, sinDV, sinV, cosV, sinB, sinD, cosD : Double;
begin
B := Latitude;
L := Longitude;
sinB := sin(B);
t := (1-e*sinB)/(1+e*sinB);
t := sqr(1+sinB)/(1-sqr(sinB)) * exp(e*ln(t));
t := k_2 * exp(alfa*ln(t));
sinU := (t-1)/(t+1);
cosU := sqrt(1-sinU*sinU);
V := alfa*L;
sinV := sin(V);
cosV := cos(V);
cosDV := cosVQ*cosV + sinVQ*sinV;
sinDV := sinVQ*cosV - cosVQ*sinV;
sinS := sinUQ*sinU + cosUQ*cosU*cosDV;
cosS := sqrt(1-sinS*sinS);
sinD := sinDV*cosU/cosS;
cosD := sqrt(1-sinD*sinD);
eps := n*arctan(sinD/cosD);
rho := rho_0*exp(-n*ln((1+sinS)/cosS));
C_X := rho*sin(eps);
C_Y := rho*cos(eps);
end;
function WGS84toJTSK(Latitude, Longitude : Double; var X, Y :double): Boolean;
// vraci JTSK souřadnice pro zadanou severní sířku a východní délku
// změřenou GPS v souřadném systému elipsoidu WGS84
var fake_alt : double;
begin
fake_alt := 0;
WGS84toBessel(Latitude, Longitude, fake_alt);
BesseltoJtsk(Latitude, Longitude, X, Y);
end;
var lat,lon,alt,x,y :double;
begin
read(lat);
read(lon);
lat := lat*pi/180;
lon := lon*pi/180;
WGS84toJTSK(lat, lon, x, y);
writeln(x:10:5, ' ', y:10:5);
end._______________________________________________
Talk-cz mailing list
Talk-cz@openstreetmap.org
http://lists.openstreetmap.org/listinfo/talk-cz
