Roman Kadaj przedstawił w 2003 podejście do rozwiązania problemu różnic pomiędzy teoretycznym a empirycznym Układem 1965, nazwane "korektami globalnymi". Jest to algorytm, dodający do wejściowych współrzędnych (teoretycznych) poprawki, aproksymowane wielomianami stopnia 16 w strefach 1-4 oraz 12 w strefie 5, których współczynniki zostały wyznaczone na podstawie różnic pomiędzy empirycznymi, a teoretycznymi współrzędnymi punktów osnowy geodezyjnej I i II klasy w każdej strefie układu. Taka korekta znacznie zbliża współrzędne teoretyczne do empirycznych, redukując różnice pomiędzy nimi zwykle do kilku centymetrów.

Wymyśliłem jak to podejście zastosować w FOSSGISie. Kod korekt globalnych Romana Kadaja napisany w Delphi przepisałem do skryptu w Pythonie. Stworzyłem w GRASSie 7.0.4 5 par map rastrowych, wyrażających w sekundach kątowych wartości przesunięć współrzędnych punktów każdej strefy PUWG 1965 pomiędzy układem latlon elipsoidy Krasowskiego, a układem latlon elipsoidy WGS84 (GRS80), dodając do tych przesunięć poprawki równe korektom globalnym policzonym skryptem corr65.py. Rozdzielczość map ustaliłem tak, by nawiązywała do zagęszczenia punktów osnowy geodezyjnej I i II klasy. Tak powstałe mapy rastrowe przekonwertowałem do gridów NTv2.  Stosując taką siatkę różnic współrzędnych zamiast obliczania przesunięcia metodą geometryczną otrzymujemy definicję dobrego przybliżenia empirycznego PUWG 1965, zrozumiałą dla PROJ.4 i oprogramowania na nim opartego. Poniżej podaję dokładny opis procesu obliczenia rozdzielczości, zasięgu przestrzennego i wartości komórek każdej mapy rastrowej oraz przetworzenia ich do postaci gridów NTv2.


1. Z serwisu geoportal.gov.pl pobrałem warstwę WFS osnowy poziomej I i II stopnia w PUWG 1992 poleceniem:

v.in.wfs wfs='http://mapy.geoportal.gov.pl/wss/service/pub/guest/G2_OSNOWA_WFS/MapServer/WFSServer?version=1.1.0&typename=G2_OSNOWA_WFS:OS_PunktOsnowyPoziomej&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=G2_OSNOWA_WFS:OS_PunktOsnowyPoziomej&SRSNAME=EPSG:2180' out=osnowa

Było to w czasach, gdy jeszcze nie wymagali autoryzacji do WFSa. Obecnie trzeba poprosić o dostęp.


2. Dla każdego punktu policzyłem odległość do jego najbliższego sąsiada:

g.copy vect=osnowa,osnowa_dist
v.db.addcol osnowa_dist col='cat_to integer','dist double'
v.distance from=osnowa_dist to=osnowa_dist from_type=point to_type=point upload=cat,dist col=cat_to,dist dmin=0.0001


3. Ponieważ parametry wielomianów użyte przez Kadaja zostały wyznaczone w oparciu o położenie punktów osnowy I i II klasy, jak sam pisze w komentarzu do swojego kodu, uznałem, że rozdzielczość gridów NTv2 powinna nawiązywać do gęstości punktów osnowy. Ustaliłem ją więc na pierwszy kwartyl odległości między dwoma najbliższymi punktami osnowy leżącymi w obrębie wszystkich arkuszy mapy topograficznej 1:10000 PUWG 1965 danej strefy. Myślę, że to sensowny kompromis między oczekiwaną dokładnością przeliczeń a szczegółowością gridów. Przyznaję, że wybrałem tą wartość „na czuja”. Jednak, jak pokazuję dalej, uzyskana dokładność przeliczeń znakomicie odpowiada danym referencyjnym, mimo bardzo niewielkich rozmiarów gridów - wszystkie 5 plików waży w sumie 424 KB. Nic nie stoi na przeszkodzie by wygenerować je w razie potrzeby od nowa używając wyższej rozdzielczości. Zasięg każdego grida wyznaczają granice arkuszy map topograficznych 1:10000 PUWG 1965, pozyskane ze skorowidza w formacie Shapefile, plus mały zapas równy rozdzielczości grida. Szczegóły:

# Set environment so that awk works properly in all languages:
unset LC_ALL; LC_NUMERIC='C'; export LC_NUMERIC
for i in 1 2 3 4 5; do
  ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:2180 -where "\"GODLO\" LIKE \"$i%.%\"" u65_10k_u92_${i}.shp 1965_B10.shp
  v.in.ogr dsn=u65_10k_u92_${i}.shp out=u65_10k_u92_${i}
  v.select atype=point btype=area alayer=1 blayer=1 ainput=osnowa_dist binput=u65_10k_u92_$i out=osnowa_dist_$i
  eval `v.univar -eg layer=1 type=point col=dist map=osnowa_dist_$i | grep '^first_quartile='`  
 
# Create a mask for each PUWG 1965 zone. The mask encompasses the extent of
# 1:10 000 PUWG 1965 zone's topo map sheets + a small buffer equal the 1st
# quartile of the distance between each 2 closest points of the geodetic control
# network within the zone's extent. This will be the target NTv2 grid extent
# per PUWG 1965 zone.
  v.buffer layer=1 type=area in=u65_10k_u92_$i out=u65_10k_u92_${i}_buf dist=$first_quartile
# Dissolve any holes within the mask, just in case.
  v.centroids in=u65_10k_u92_${i}_buf out=u65_10k_u92_${i}_buf_addcat opt=add
  v.dissolve in=u65_10k_u92_${i}_buf_addcat out=u65_10k_u92_${i}_buf_diss layer=1
# Create a buffer of the central point of the region encompassing the mask
# created in previous step. The buffer radius is half the $first_quartile. This
# will be the basis for calculating the NTv2 grid's resolution in EPSG:4179
# location.
  g.region vect=u65_10k_u92_${i}_buf
  v.in.region out=u65_10k_u92_${i}_buf_reg
  v.buffer type=centroid in=u65_10k_u92_${i}_buf_reg out=u65_10k_u92_${i}_buf_reg_res dist=`echo $first_quartile | awk '{printf "%.15f", $1/2.0}'` --o
done

W kodzie powyżej wyznaczam rozdzielczość grida poprzez stworzenie mapy wektorowej składającej się tylko z bufora wokół centralnego punktu grida w układzie metrycznym i transformacji tego obiektu do układu latlon, ponieważ w układzie latlon długość osi poziomej jest zmienna w zależności od położenia punktu na osi pionowej – zmniejsza się w kierunku biegunów, a zwiększa w kierunku równika. Z tego powodu ustalenie rozdzielczości grida w układzie latlon jest nieco utrudnione. Łatwiej mi było to policzyć w „symetrycznym” układzie metrycznym i tak uzyskany wynik transformować do układu latlon.


4. Mapy wektorowe zasięgów gridów w PUWG 1992 konwertuję do układu latlon na elipsoidzie Krasowskiego (EPSG 4179):

for i in 1 2 3 4 5; do
  v.proj in=u65_10k_u92_${i}_buf_diss location=gsb_u92 mapset=wyk1
  v.proj in=u65_10k_u92_${i}_buf_reg_res location=gsb_u92 mapset=wyk1
done


5. W lokacji EPSG 4179 rasteryzuję mapy wektorowe wyznaczające zasięg każdego grida w specyficznej dla niego rozdzielczości, która jest pochodną gęstości punktów osnowy geodezyjnej I i II stopnia w danej strefie Układu 65 (jak opisałem w #3):

# Set environment so that awk works properly in all languages:
unset LC_ALL; LC_NUMERIC='C'; export LC_NUMERIC
# Rasterize each vector "mask" in a given zone's *resolution* and extent:
for i in 1 2 3 4 5; do
  eval `v.info -g map=u65_10k_u92_${i}_buf_reg_res`
  nsres=`echo $north $south | awk '{printf "%.15f", $1-$2}'`
  ewres=`echo $east $west | awk '{printf "%.15f", $1-$2}'`
  for j in 1 2 3 4 5; do
    g.region -a vect=u65_10k_u92_${j}_buf_diss nsres=$nsres ewres=$ewres
    v.to.rast type=area layer=1 in=u65_10k_u92_${j}_buf_diss out=u65_10k_u92_${j}_buf_diss_${i} use=val val=1
  done
done


6. Ostatni i najbardziej złożony etap – wytłumaczenie w komentarzach wewnątrz kodu.

Proszę zwrócić uwagę, że zamiast domyślnych 7 parametrów przesunięcia między datum Pułkowo elipsoidy Krasowskiego a WGS84, tak jak zdefiniowano je w bazie danych EPSG, używam wartości o większej precyzji. Pochodzą one z Instrukcji Technicznej G-2 (wyd 5. zmienione, GUGiK 2001, Warszawa) s. 29. Takie same wartości podawał Roman Kadaj w znanych mi publikacjach. Próba zastosowania mniej precyzyjnych, domyślnych parametrów EPSG do transformacji współrzędnych między elipsoidą Krasowskiego a WGS84 powoduje systematyczne przesunięcie około 2 cm względem wyników referencyjnych podanych w pliku przyklady.txt, publikowanym kilka lat temu na geonet.net.pl wraz z plikiem kod_korekta65.txt zawierającym kod korekt globalnych. (Ze względu na usunięcie plików z ich oryginalnej lokalizacji udostępniam je w moim repozytorium w celach porównawczych i archiwalnych. Nie roszczę sobie żadnych praw do zawartości tych plików). Taki błąd można zignorować np. dla skanów map topograficznych 1:10000, gdzie pojedynczy piksel może mieć kilka metrów bez straty dla czytelności i rzetelności mapy, ale w przypadku korekt globalnych walczymy o centymetry, więc użyłem dokładniejszych parametrów, podanych przez GUGiK.

Dla pełniejszego oglądu sprawy dodam, że jakieś 10 lat temu, kiedy korespondowałem z EPSG w celu poprawienia i uzupełnienia kilku definicji polskich układów w ich bazie danych, próbowałem przekonać ich, by zamiast tych zaokrąglonych w dół wartości parametrów towgs84 pozyskanych od Niemców użyli ich dokładniejszych odpowiedników z G-2. Odmówiono mi, argumentując jak następuje: „I (…) found the difference to be 2 cm +/– 5 mm. In itself it is questionable whether this is sufficient to make it worthwhile changing the database. (…) Instead of introducing a new record I would suggest adding some remarks to our existing record. But there is another technical reason (...). At this greater precision the assumption that the method can be used for the reverse transformation by changing the signs of the parameter values is no longer valid, necessitating separate transformations for Pulkovo 1942(58) > ETRS89 and ETRS89 > Pulkovo 1942(58)” oraz „The tests by EuroGeographics suggest that the transformation for Poland has, on average, an accuracy of about 1m. Adding extra resolution (more digits) to the transformation parameter values will not increase the accuracy of applying the transformation”.

Mimo to ewidentny pozostaje fakt, że tylko zastosowanie dokładniejszych parametrów podanych przez GUGiK pozwala wygenerować gridy realizujące dwustronną transformację Krasowski ← → WGS84 z różnicami względem danych referencyjnych Kadaja maksymalnie 2 mm (patrz #8 poniżej). Użycie zgeneralizowanych współczynników towgs84, jakie obecnie są przypisane w bazie EPSG do polskich układów opartych na elipsoidzie Krasowskiego, skutkuje powiększeniem tego błędu o około 2 cm. Za skorzystaniem z parametrów podanych w G-2 przemawia także zgodność wyników uzyskanych przy ich użyciu z narzędziami przeliczeń współrzędnych GDAL i PROJ.4 z wynikami transformacji teoretycznej uzyskanymi w TRANSPOLU 2.06. Jak należało się spodziewać różnią się one od analogicznych wyników dla zaokrąglonych w dół współczynników od EPSG za EuroGeographics o rzeczone 2 cm - patrz Tab. 1.




















Tab. 1. Porównanie wyników przeliczeń punktów referencyjnych BLH GRS80 z przyklady.txt za pomocą transformacji teoretycznej w TRANSPOLu 2.06 z wynikami uzyskanymi cs2cs z użyciem domyślnej definicji towgs84 EPSG dla polskich układów na elipsoidzie Krasowskiego oraz definicji towsg84 o wyższej precyzji, pochodzącej z Instrukcji Technicznej G-2.


Podsumowując: Gdyby ktoś chciał przeliczać za pomocą moich gridów między PUWG 1965 a innym polskim układem opartym na elipsoidzie Krasowskiego (PUWG 1942, GUGiK 80, latlon/Krasowski – EPSG:4179), albo między „teoretycznym” a „empirycznym” Układem 1965 proszę pamiętać, że definicje EPSG tych układów, dostarczane np. razem z bibliotekami GDAL, PROJ.4, GeoTIFF, zawierają mniej precyzyjne 7 parametrów towgs84 od EuroGeographics. Zatem aby uniknąć 2 cm przesunięcia należy je zastąpić tymi z G-2, co z kolei uniemożliwia posługiwanie się wprost kodami EPSG tych układów. Prostym rozwiązaniem problemu jest rozwinięcie kodu EPSG do jego jawnej definicji w formacie PROJ.4 (np. za pomocą epsg_tr.py lub gdalsrsinfo) i ręczna modyfikacja parametrów towgs84. Jak np. robię to w kodzie poniżej.

# Set environment so that awk works properly in all languages:
unset LC_ALL; LC_NUMERIC='C'; export LC_NUMERIC
# Create a global mask at each 1-5 zone's specific resolution, and trim each
# zone's grid to the mask:
_mkgsb() {
  # This function calls corr65.py (Python implemenentation of Kadaj's global
  # corrections algorithm) to translate the location of cells within each PUWG
  # 1965 zone grid mask, expressed in Krassovsky latlon coordinates converted to
  # PUWG 1965, into WGS84 latlon coordinates + the global corrections.
  # Substracting the two gives the final product - NTv2 transformation grids,
  # which have the Kadaj's transformation coefficients "embedded" in.
  nam=$1; i=$2
  # Instead of EPSG codes for PUWG 1965 zones I'm using explicit PROJ.4
  # definitions in order to override the default EPSG's
  # towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 with the values officially
  # recommended by the Polish state geodesy authority (GUGiK) in their "G-2"
  # guidelines. It might be argued that this doesn't necessarily improve real-
  # life accuracy of the transformation. Anyway, I'm doing this to avoid ~2 cm
  # systematic difference between PUWG 1965 coordinates obtained with default EPSG
  # towgs84 coefficients and the reference, benchmark PUWG 1965 coordinates
  # published in przyklady.txt on geonet.net.pl, which apparently were obtained
  # using those official, higher-precision coefficients.
  if [ $i -eq 1 ]; then epsg_to='+proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +no_defs'; fi
  if [ $i -eq 2 ]; then epsg_to='+proj=sterea +lat_0=53.00194444444445 +lon_0=21.50277777777778 +k=0.9998 +x_0=4603000 +y_0=5806000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +no_defs'; fi
  if [ $i -eq 3 ]; then epsg_to='+proj=sterea +lat_0=53.58333333333334 +lon_0=17.00833333333333 +k=0.9998 +x_0=3501000 +y_0=5999000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +no_defs'; fi
  if [ $i -eq 4 ]; then epsg_to='+proj=sterea +lat_0=51.67083333333333 +lon_0=16.67222222222222 +k=0.9998 +x_0=3703000 +y_0=5627000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +no_defs'; fi
  if [ $i -eq 5 ]; then epsg_to='+proj=tmerc +lat_0=0 +lon_0=18.95833333333333 +k=0.999983 +x_0=237000 +y_0=-4700000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +no_defs'; fi
  espg4179='+proj=longlat +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +no_defs'
  # Create raster maps of Krassovsky longituteds and lattitudes in a PUWG 1965
  # zone:
  r.mapcalc expression="lon_frK_${nam}=double(x())"
  r.mapcalc expression="lat_frK_${nam}=double(y())"
  # Use the latlon coorinates stored in those maps as (indirect) input for
  # Kadaj's algorithm:
  r.stats -1n in=lon_frK_$nam,lat_frK_$nam | while read lonlat_frK; do
    echo $lonlat_frK | gdaltransform -s_srs "$espg4179" -t_srs "$epsg_to" | awk '{print $2,$1}' | while read yx_K; do
      lonlat_to=`corr65.py $yx_K $i E | awk '{print $2,$1}' | gdaltransform -s_srs "$epsg_to" -t_srs EPSG:4326 | awk '{print $1,$2}'`
      echo $lonlat_frK $lonlat_to
    done
  done > frK_to.txt
  r.in.xyz --o in=frK_to.txt out="lon_to_$nam" z=3 fs=space type=DCELL
  r.in.xyz --o in=frK_to.txt out="lat_to_$nam" z=4 fs=space type=DCELL
  # K+x=G, x=G-K:
  # https://github.com/Esri/ntv2-file-routines/blob/master/README.md:
  #
  # "Longitudes are reversed from normal usage, with values
  # west of Greenwich being positive, and values east of Greenwich
  # being negative. This is probably because all of Canada (the
  # creators of this file format) is west of Greenwich":
  r.mapcalc expression="lat_diff_${nam}=double(lat_to_${nam}-lat_frK_${nam})*double(3600.0)"
  r.mapcalc expression="lon_diff_${nam}=-double(lon_to_${nam}-lon_frK_${nam})*double(3600.0)"
  # Replace nulls with 0s. Otherwise the features overlapping the null areas
  # would get lost in translation:
  r.null map=lat_diff_$nam null=0
  r.null map=lon_diff_$nam null=0
  # And remove the mask, or the 0s will be exported as nulls regardless:
  g.remove type=rast name=MASK -f
  # Export the difference = the NTv2 GSB grid. NOTE: r.out.ntv2 does these 2
  # steps in one batch:
  i.group group=diff_$nam in=lat_diff_$nam,lon_diff_$nam
  # Notice we are converting from double to floating point precision - NTv2
  # requirement.
  r.out.gdal -f input=diff_$nam format=NTv2 out=${nam}.gsb type=Float32
}

for i in 1 2 3 4 5; do
  eval `v.info -g map=u65_10k_u92_${i}_buf_reg_res`
  nsres=`echo $north $south | awk '{printf "%.15f", $1-$2}'`
  ewres=`echo $east $west | awk '{printf "%.15f", $1-$2}'`
  # Create a mask at the given zone's own resolution, but with a whole PUWG 1965
  # extent:
  g.remove type=rast name=MASK -f
  g.region -a rast=u65_10k_u92_1_buf_diss_${i},u65_10k_u92_2_buf_diss_${i},u65_10k_u92_3_buf_diss_${i},u65_10k_u92_4_buf_diss_${i},u65_10k_u92_5_buf_diss_${i} nsres=$nsres ewres=$ewres
  r.patch in=u65_10k_u92_1_buf_diss_${i},u65_10k_u92_2_buf_diss_${i},u65_10k_u92_3_buf_diss_${i},u65_10k_u92_4_buf_diss_${i},u65_10k_u92_5_buf_diss_${i} out=u65_10k_u92_all_buf_diss_${i}
  g.region -a nsres=$nsres ewres=$ewres zoom=u65_10k_u92_${i}_buf_diss_${i}
  r.mask rast=u65_10k_u92_all_buf_diss_${i}
  _mkgsb u65_10k_u92_${i}_buf_diss $i
done

Tak powstałe 5 gridów wrzuciłem do repo. Przy okazji napisałem skrypt r.out.ntv2 w Pythonie ułatwiający eksport gridów NTv2 z GRASSa. Udostępniam go w repozytorium GRASS 7 addons.


7. Teraz zobaczmy jak gridy spisują się w praktyce. Format NTv2 jest zrozumiały dla PROJ.4, a przez to dla GDAL/OGR, GRASS, QGIS itp. Aby użyć grida transformacyjnego należy zastąpić 7-parametrową definicję towgs84 definicją nadgrids, podając w niej ścieżkę do pliku w formacie NTv2.

7.1. Weźmy dla przykładu konwersję współrzędnych punktu w okolicach Lublina między WGS84 (EPSG:4326) a PUWG 1965 (Lublin leży w strefie I – więc EPSG:3120).

Sprawdźmy jak wyglądają definicje obu układów w formacie PROJ.4:

$ gdalsrsinfo -o proj4 EPSG:3120
'+proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs '
$ gdalsrsinfo -o proj4 EPSG:4326
'+proj=longlat +datum=WGS84 +no_defs'
„datum=WGS84” znaczy tyle, co „+ellps=WGS84 +towgs84=0,0,0”. Zatem kodu EPSG:4326 możemy używać wprost. Natomiast w przypadku EPSG:3120 należy podmienić definicję przesunięcia towgs84 na bardziej szczegółową z Instrukcji Technicznej G-2 aby uniknąć omawianego wcześniej systematycznego błędu ~2 cm. No to przeliczamy. Posłużę się narzędziem cs2cs z PROJ.4 (gdaltransform z GDAL dałby równie dobre wyniki).
$ echo '22.57 51.24' | cs2cs -f '%.5f' +init=epsg:4326 +to +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m
4740923.25402 5536481.18439 -28.89913

I w drugą stronę dla pewności, że transformacja jest odwracalna:

$ echo '4740923.25402 5536481.18439 -28.89913' | cs2cs -f '%.5f' +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +to +init=epsg:4326
22.57000 51.24000 0.00009

7.2. Teraz wykonajmy te same przekształcenia zastąpiwszy parametr towgs84 przez nadgrids ze ścieżką do stworzonego jak wcześniej opisałem grida NTv2, uzyskując w ten sposób definicję empirycznego PUWG 1965, tj. zawierającą korekty globalne a'la Kadaj:

$ echo '22.57 51.24' | cs2cs -f '%.5f' +init=epsg:4326 +to +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +nadgrids=~/u65_10k_u92_1_buf_diss.gsb +units=m
4740923.29650 5536481.05934 0.00000

Tu też sprawdzamy odwracalność transformacji:

$ echo '4740923.29650 5536481.05934 0.00000' | cs2cs -f '%.5f' +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +nadgrids=~/u65_10k_u92_1_buf_diss.gsb +units=m +to +init=epsg:4326
22.57000 51.24000 0.00000

Uzyskana różnica około 4 cm w poziomie i 12 cm w pionie między współrzędnymi PUWG 1965 w #7.1 („teoretycznymi”) oraz 7.2 („empirycznymi”) to właśnie efekt zastosowania korekt globalnych.

7.3. Dla porównania przeliczmy wzajemnie teoretyczne współrzędne płaskie PUWG 1965 uzyskane w #7.1 i empiryczne z #7.2, używając grida:

$ echo '4740923.25402 5536481.18439' | cs2cs -f '%.5f' +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m +to +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +nadgrids=~/u65_10k_u92_1_buf_diss.gsb +units=m
4740923.29716 5536481.05941 28.89920
$ echo '4740923.29650 5536481.05934' | cs2cs -f '%.5f' +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +nadgrids=~/u65_10k_u92_1_buf_diss.gsb +units=m +to +proj=sterea +lat_0=50.625 +lon_0=21.08333333333333 +k=0.9998 +x_0=4637000 +y_0=5467000 +ellps=krass +towgs84=33.4297,-146.5746,-76.2865,-0.35867,-0.05283,0.84354,-0.84077 +units=m
4740923.25403 5536481.18439 -28.89913

Na koniec te same współrzędne przeliczmy z teoretycznego do empirycznego PUWG 1965 i z powrotem skryptem corr65.py, implementującym korekty globalne:

Usage:
corr65.py x y n r

Where:
x: source PUWG 1965 x coordinate
y: source PUWG 1965 y coordinate
n: PUWG 1965 zone number (1, 2, 3, 4 or 5)
r: source coordinates type - "T"heoretical, "E"mpirical

Example:
Converting theoretical coordinates to empirical in zone 1:
corr65.py 5467000 4637000 1 T

Z teoretycznego do empirycznego:

$ corr65.py 5536481.18439 4740923.25402 1 T
5536481.05926 4740923.2971

I w odwrotną stronę:

$ corr65.py 5536481.05934 4740923.29650 1 E
5536481.18449 4740923.25353

7.4. Jak widać w powyższych przykładach zachodzi duża zgodność przeliczeń współrzędnych płaskich wykonanych skryptem corr65.py, realizującego algorytm korekt globalnych Kadaja, i cs2cs z gridem NTv2, wygenerowanym z pomocą tegoż skryptu. Zwraca uwagę około 28 m przesunięcie wysokości w przeliczeniach cs2cs z gridem z #7.3. Wynosi ono tyle samo co przy przeliczaniu między elipsoidą Krasowskiego a WGS84 w #7.1. O ile w tamtym przypadku różnica wysokości jest spodziewana ze względu na odstępy pomiędzy elipsoidami różnych rozmiarów (porównaj np. „Elipsoidy a układy” cz. 5 na geoforum.pl), to w przypadku transformacji pomiędzy dwoma układami opartymi na tej samej elipsoidzie może dziwić. Tłumaczę to tym, że 7-parametrowa transformacja towsgs84 modyfikuje współrzędne w trójwymiarze, podczas gdy grid NTv2 dostarcza jedynie wartości przesunięć na płaszczyźnie. Tym samym przesunięcia w osiach X i Y aplikowane przez transformacje towsgs84 i nadgrids wzajemnie się znoszą, pozostawiając jedynie przesunięcie o wartość równą korektom globalnym Kadaja, podczas gdy przesunięcie wysokości pozostaje. Jednak rzutuje ono na przeliczenia współrzędnych płaskich w sposób zazwyczaj pomijalny, generując odchyłki w okolicach zaledwie 1 mm, jak obserwujemy również w tym przypadku.

Narzędziom GDAL/OGR jak gdalwarp, ogr2ogr należy podawać układ współrzędnych z gridem NTv2 w takim sam sposób jak dla cs2cs.

W QGISie podobnie. Należy zdefiniować 5 układów współrzędnych użytkownika w formacie PROJ.4, o parametrach jak w EPSG:3120, 2172-2175, ale z sekcją towgs84 zastąpioną przez nadgrids analogicznie do powyższych przykładów.

W GRASSie edytujemy plik PROJ_INFO w mapsecie PERMANENT danej lokacji, zastępując linię „towgs84: 33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84” linią „nadgrids: nazwa_pliku_grida_NTv2”. Plik grida musi zostać ręcznie skopiowany do katalogu $GISBASE/etc/proj/nad/. Tak jest w GRASSie 6.x i 7.0.x. Od wersji 7.x prawdopodobnie grid będzie musiał znaleźć się w odpowiednim katalogu PROJ.4 - najczęściej /usr/share/proj/.


8. Jak widać w Tab. 2 porównanie wyników obliczeń moją implementacją korekt globalnych w postaci gridów NTv2 z danymi referencyjnymi w pliku przyklady.txt wypada pozytywnie. W skrajnych przypadkach różnica współrzędnych sięga zaledwie 2 mm. Proszę zignorować różnicę kilkudziesięciu metrów dla punktu nr 40 w strefie IV i V. Zasięg moich gridów wykracza tylko nieco poza prostokąt obejmujący arkusze mapy topograficznej 1:10000 danej strefy PUWG 1965, a punkt nr 40 nie leży w obrębie tych zasięgów dla strefy IV i V. Przy okazji widzimy co się dzieje, gdy próbujemy za pomocą grida NTv2 transformować współrzędne punktu nie pokrywającego się z jego zasięgiem – PROJ.4 po cichu aplikuje zerowe przesunięcie. Oczywiście nie ma technicznych przeszkód by wygenerować gridy dla dowolnie dużego obszaru. Ja jednak przyjąłem, jak wydaje mi się – sensowne, ograniczenie do zasięgu arkuszy 1:10000 z małym marginesem i na razie tego się trzymam.

przyklady3.png




























Rys 1. Rozmieszczenie punktów kontrolnych z przyklady.txt na tle orientacyjnych granic 5 gridów NTv2 i zasięgu arkuszy mapy 1:10000 PUWG 1965.




























Tab. 2. Porównanie referencyjnych współrzędnych empirycznych PUWG 1965 z przyklady.txt z wynikami przeliczeń tychże punktów z BLH GRS80 na PUWG 1965 z użyciem moich gridów NTv2.



Trzeba pamiętać, że korekty globalne zazwyczaj nie są panaceum na zniekształcenia PUWG 1965. Zarówno Kadaj jak i inni autorzy (np. http://www.infraeco.pl/pl/art/a_16170.htm?plik=932, http://www.numerus.net.pl/transformacja.html) wskazują na jego ograniczoną przydatność w niektórych zastosowaniach i częstą konieczność dodatkowego użycia korekt lokalnych, jeżeli zależy nam na dużej dokładności.

Myślę jednak, że dobrze mieć korekty globalne PUWG 1965 w arsenale FOSSGIS. Wiele krajowych narzędzi już dawno zaimplementowało to rozwiązanie - np. SWDE-konwerter, GEONET, EWMAPA, GEO-TRANS, GEO-INFO, GeoKonwerter, C-GEO. Jak rozumiem nie ma po temu formalnych przeszkód. Roman Kadaj napisał w komentarzu do swojego kodu korekt globalnych, niegdyś do pobrania ze strony geonet.net.pl, co następuje: "W celu zachowania jednolitości opracowań geodezyjnych i kartograficznych (dotyczy to zwłaszcza systemów EGiB), z inicjatywy Departamentu Katastru GUGiK, dokonuje się upowszechnienia tej procedury i jej parametrów jako ogólnodostępnej informacji technicznej (...). Niniejsza publikacja autorska służy również temu celowi. Przy wszelkich zastosowaniach należy zwrócić szczególną uwagę na jej bezbłędne skopiowanie", a we wpisie na geonet.net.pl 24.03.2003 p.t. "Korekty "1965" opublikowane" napisano "Z chwilą publikacji nie ogranicza się jakichkolwiek zastosowań".

Kiedy kończyłem pracę nad moją implementacją korekt globalnych Układu 1965 natknąłem się na publikację Tomasza Świętonia, zamieszczoną na geonet.net.pl p.t. „Ocena możliwości zastosowania siatek korekt lokalnych do opracowania jednolitego algorytmu transformacji pomiędzy układami 1965 a 1992 na terenie Polski”. Ku mojej radości dowiedziałem się od autora, że w rezultacie dalszych prac nad tym zagadnieniem TRANSPOL 2.06 implementuje już od dłuższego czasu korekty empiryczne za pomocą gridów. I to nie tylko korekty globalne dla samego Układu 1965, które są zaledwie wstępem do uzyskania dokładności geodezyjnej i ograniczone do tego jedynie układu, ale także bardziej zaawansowane korekty empiryczne wysokiej precyzji dla wszystkich polskich układów współrzędnych. Umożliwia to wykorzystanie wyników przeliczeń uzyskanych TRANSPOLem w geodezyjnych opracowaniach technicznych. Idzie nie tylko o korekty dla układów opartych na elipsoidzie Krasowskiego. Także PL-ETRF89 i PL-ETRF2000 na GRS80, będące podstawami do realizacji układów współrzędnych płaskich 1992, 2000 i UTM w różnych okresach, charakteryzują się niesystematycznymi różnicami rzędu kilku cm, wymagającymi zastosowania odpowiednich poprawek. TRANSPOL nie używa jednak gridów NTv2, a stosuje połączenie kilku gridów we własnościowym formacie ze specyficznymi dla różnych układów algorytmami przeliczeń na bazie tych gridów. Mimo to powinno dać się go użyć jako źródła współrzędnych potrzebnych do obliczenia przesunięć między wybranymi parami układów w podobny sposób, jak użyłem corr65.py w kodzie cytowanym powyżej. Pozwala to na stworzenie zestawu gridów NTv2 realizujących te same przeliczenia, co TRANSPOL 2.06 z włączoną opcją korekt empirycznych. Kiedyś może spróbuję ugryźć ten temat. Takie gridy mogą ułatwić wykorzystanie GRASSa, GDALa i QGISa w wymagających wysokiej precyzji pracach geodezyjnych na terenie Polski.