SWDB - dane wodne jako dopasowanie SRTM

[Ra]

Aby wypozycjonować (skalibrować) dane SRTM dla mapy topograficznej z siatką kilometrową, postanowiłem zacząć od danych wodnych SWDB (kształty zbiorników wodnych), dołączanych do SRTM. (Rzecz dotyczy aktualnie scenerii PMP-PW, ale może się to przydać dla innych rejonów też.)

Napisałem sobie proste klasy w C++, które wczytują tzw. shapefile (.shp), w którym zapisane są wektorowo kształty zbiorników. (Dodatkowo są dwa pliki - indeksowy oraz .DBF - na razie nieprzydatne.) Shapefile ma troszkę zamotaną strukturę, z uwagi na pomieszane Big Endian z Little Endian oraz długości podawane w słowach (dwubajtach).

Wektorowe kształty zbiorników to rażąca pikseloza. Wygląda to tak, jakby ktoś na zdjęciu satelitarnym o niskiej rozdzielczości zaznaczył "niebieskie" obszary narzędziem wypełniającym. Piksele te zostają obrysowane łamaną o kątach 90° - dane zupełnie nieprzydatne do dalszego przetwarzania, chyba żeby aproksymować łamaną za pomocą krzywej Béziera. Współrzędne punktów są zapisane w stopniach.

Aby narysować te obrysy na mapie, użyłem funkcji przekształcenia:

x_m=a*x+b*y+c*x*y+d
y_m=e*x+f*y+g*x*y+h

Następnie dopasowywałem współczynniki tak, aby obrysy pokryły się ze zbiornikami Pogoria III, Zalew Sosina, J. Dziećkowice. Udało się to osiągnąć przy następujących współczynnikach:
a=71350 - [m/°] w poziomie
b=2850 - obrót
c=-0.6544 -skalowanie
d=1510770 - przesunięcie
e=1600 - obrót
f=-111200 - [m/°] w pionie
g=-0.05179 -skalowanie
h=5556800 - przesunięcie

Z kolei dane wyliczone z map topograficznych wyglądają następująco:
- odległość między 19°E a 20°E wzdłuż równoleżnika 51°N: 70406m,
- odległość między 19°E a 20°E wzdłuż równoleżnika 50°N: 71695m,
- odległość między 50°N a 51°N wzdłuż południka 19°E: 111339m,
- odległość między 50°N a 51°N wzdłuż południka 20°E: 111237m.

W efekcie wspomniane trzy zbiorniki pasowały prawie idealnie. Zbiorniki wodne w okolicy Sosnowiec-Mysłowice pasowały w miarę dobrze. Dla porównania zrobiłem sobie mapkę okolicy Koniecpola (13 km od Włoszczowy w stronę Częstochowy), z zawartym tam zbiornikiem. Zbiornik ten odstawał ok. 4km na północny-wschód. Wszelkie próby sprowadzenia go na wskazane na mapie miejsce kończyły się znacznymi odchyleniami wcześniej ustawionych trzech zbiorników.

Dalsze możliwości są następujące:
- ograniczyć przetwarzanie SRTM do obszaru, w którym dopasowanie danych wodnych ma dosyć dużą dokładność (nie chodzi tylko o SRTM, bo można wykorzystać inne dane, tworzone w oparciu o GPS - stosując tę samą konwersję współrzędnych, co dla SWDB);
- wybrać punkty charakterystyczne z SWDB (np. końce zbiorników, ostre łuki) i wyliczyć ich docelową pozycję na mapie; następnie wyliczać współczynniki metodą najmniejszych kwadratów/sześcianów;
- zmienić funkcję konwersji - np. przeliczanie współrzędnych biegunowych.

Pokażę mniej więcej, jak to wygląda. Czerwonym kolorem są tory (wstępnie przesunięte, ale nie dopieszczone - widać ząbek po lewej na dole), a jasnoniebieskim jest zarys zbiorników z SWDB (załącznik 1).

Dla porównania, okolice Koniecpola z zaznaczonym zbiornikiem (załącznik 2). Jest też niewielkie prawdopodobieństwo, że ten zbiornik jest błędnie narysowany... trzeba by znaleźć kolejny i próbować dalej.

Znalazłem jakieś wzory przekształcające współrzędne geograficzne na współrzędne prostokątne płaskie. Nie wyglądają ciekawie... http://www.navi.pl/gps/artykuly/uklady/uklady.html.

[Ra]

No i niestety wprowadziłem się w błąd. Ten duży zbiornik koło Koniecpola nie jest tym, za który go uważałem. Na zrzucie ekranu okolicy Koniecpola, na dole, jest widoczny mały zbiornik - to on powinien być przesunięty o kilkaset metrów w lewo (na zachód). Tak wynika z przeliczenia współrzędnych, które znalazłem.

Dostałem odpowiedź na forum dla geodetów, abym poszukał na stronie Tadeusza Syryjczyka - http://www.syryjczyk.krakow.pl/Mapy_Polskie_GPS.htm. Czytałem kiedyś zawartość tej strony, ale wtedy nie zwróciłem uwagi na zamieszczony tam kalkulator Edwarda Zatorskiegow postaci arkusza XLS.

Interesuje mnie przeliczenie ze współrzędnych geograficznych WGS84 na współrzędne prostokątne płaskie Pułkowo 1942, więc do tego się ograniczę, chociaż kalkulator umożliwia również inne przeliczenia.

Obliczenia wykonywane są na trzeciej zakładce arkusza. Dane w stopniach minutach i sekundach wpisane są na drugiej zakładce, w pola D11, F11, H11 (B-szerokość) oraz D12, F12, H12 (L-długość). Na trzeciej zakłace są przeliczane na stopnie w polu C135 (B) oraz C136 (L).

Wynik w postaci punktu (X,Y) [m] we współrzędnych prostokątnych płaskich wyliczany jest do pól C142 (X) oraz C143 (Y). Dla interesującego mnie przekszałcenia (B148=10, B149=6, B152=2), dane te pobierane są z pól H117 (X) oraz H118 (Y). Tam są następujące formuły:

H117=H77*H104-0
H118=GDY(B152=1;H77*H105+3500000;H77*H105+4500000)

Z kolei:

H77=1

H104=H86*(H96+(H98*SIN(2*H96)*COSH(2*H97))
+(H99*SIN(4*H96)*COSH(4*H97))
+(H100*SIN(6*H96)*COSH(6*H97))
+(H101*SIN(8*H96)*COSH(8*H97)))

H105=H86*(H97+(H98*COS(2*H96)*SINH(2*H97))
+(H99*COS(4*H96)*SINH(4*H97))
+(H100*COS(6*H96)*SINH(6*H97))
+(H101*COS(8*H96)*SINH(8*H97)))

H86=6367558.49687
H97=0.5*LN((1+COS(H80)*SIN(H81-H84))/(1-COS(H80)*SIN(H81-H84)))
H98=0.0008376117571403
H96=ATAN(SIN(H80)/(COS(H80)*COS(H81-H84)))
H99=0.0000007606346141534
H100=0.000000001197122824063
H101=0.000000000002441972616146

H80=H93=2*(ATAN(POTĘGA((1-H88*SIN(H82))/(1+H88*SIN(H82));H88/2)*TAN(H82/2+H85/4))-H85/4)
H81=RADIANY(H79)
H84=RADIANY(H76)

H88=0.0818133340169
H82=RADIANY(H73)
H85=PI
H79=H74=D58=GDY(B148<9;C58;GDY(B148=9;C136;L61))
H76=GDY(B152=1;15;GDY(B152=2;21;24)) - południk odniesienia

H73=D57=L60

C57=0
C58=0

L60=STOPNIE(L59)
L61=STOPNIE((ACOS(L38/L49)))
L59=ATAN((L40+L58)/L49)
L38=L34=L19+L22*L19+L23*L20+L24*L21+L31
L49=POTĘGA((POTĘGA(L38;2)+POTĘGA(L39;2));1/2)
L40=L36=L17+L28*L19+L29*L20+L30*L17+L33
L58=L42*L45*L57/POTĘGA(1-L57*L57;1/2)
L39=L35=L20+L25*L19+L26*L20+L27*L21+L32
L19=L15=(L10+L7)*COS(L8)*COS(L9)
L22=0.00000084007644
L23=0.00000408960694
L20=L16=(L10+L7)*COS(L8)*SIN(L9)
L24=0.00000025613907
L21=L17=(L10*(1-L14)+L7)*SIN(L8)
L31=-33.4297
L42=6378245
L45=0.0818133340169
L57=L45*SIN(L56)
L28=-0.00000025614618
L29=0.00000173888682
L30=0.00000084077125
L33=76.2865
L25=-0.0000040896065
L26=0.00000084076292
L27=-0.00000173888787
L32=146.5746
L10=L11/POTĘGA(1-L14*SIN(L8)*SIN(L8);1/2)
L7=200
L8=RADIANY(L5)
L9=RADIANY(L6)
L56=ATAN((L40+L55)/L49)
L14=L13*L13
L11=6378137
L55=L42*L45*L54/POTĘGA(1-L54*L54;1/2)
L5=C135
L6=C136
L13=0.0818191910428
L54=L45*SIN(L53)
L53=ATAN((L40+L52)/L49)
L52=L42*L45*L51/POTĘGA(1-L51*L51;1/2)
L51=L45*SIN(L50)
L50=ATAN((L40+0)/L49)

Formuły te wpisałem do programu w C++, zmieniając ATAN(x/y) na ATAN2(x,y), POTĘGA(x,1/2) na sqrt(x) itd.

Wyniki obliczeń wychodzą mi identyczne, jak te w arkuszu. Dane wodne, przekształcone tymi wzorami wyglądają tak, jakby powinny być przesunięte o 10m w lewo i 10m w górę - ale to może już wynikać z niedokładności ich pomiaru. W każdym bądź razie, moje wcześniejsze przekształcenie było całkiem dobre.