3  Opšte jednačine kartografskih projekcija

U narednom poglavlju predstavljeni su ključni koncepti kartografskog preslikavanja, praktično ovi koncepti se mogu primeniti na bilo koju konkretnu projekciju. Ovo poglavlje obuhvata sledeće celine:

• Razmatranje parametarskih jednačina sfere i parametarskih jednačina elipsoida kako bismo bolje razumeli osnovne geometrijske modele kod funkcija kartografskog preslikavanja.

• Gausove fundamentalne veličine koje su ključne za razumevanje neizbežnih deformacija koje prate kartografske projekcije prilikom preslikavanja krivolinijskih koordina u pravougle ravanske koordinate.

• Proučavanje opšte teorije kartografskog preslikavanja kako bismo saznali kako se prostorni podaci transformišu sa trodimenzionalnog prostora na dvodimenzionalne karte.

• Detaljan opis linijskog elementa u ravni karte, linearne razmere, razmere površina i deformacija uglova.

• Koncept i računanja koja se odnose na elipsu deformacija, Tisoovu indikatrisu, pružajući uvid u karakteristike deformacija pri preslikavanju sa sfere ili elipsoida na ravne karte.

• Podela projekcija prema karakteru deformacija, klasifikacija projekcija u odnosu na pomoćne projekcione površine i položaj pomoćnih projekcionih površina, kao i klasifikacija projekcija u odnosu na zone.

Kroz ove celine čitaoci će razviti dublje razumevanje osnovnih koncepata kartografskog preslikavanja, što će im omogućiti da se upuste u izvođenje formula i razumevanje bilo koje konkretne kartografske projekcije koja se koristi u radu sa prostornim podacima.

3.1 Parametarske jednačine sfere

U prethodnom poglavlju smo diskutovali o parametrizaciji sfere, ovde ćemo još jednom ponoviti. Meridijani i paralele na sferi predstavljaju parametarske, tj. koordinatne linije, to je jedna od mogućih parametrizacija sfere. Pozicija bilo koje tačke na sferi može se definisati presekom jednog meridijana i jedne paralele. Kod slučaja kada smo imali ravan gde je bilo koja tačka određena presekom samo jedne x i samo jedne y koordinatne linije, to je slučaj regularne parametrizacije površi. Kod sfere imamo slučaj da postoje i singularne tačke, a to su polovi u kojima se sustiče svaki meridijan, a poluprečnik paralele u polovima je 0. Polovi nisu određeni presekom samo jedne paralele i samo jednog meridija, zato su to singularne tačke, pa i parametrizacija sfere meridijanama i paralelama nije u potpunosti regularna parametrizacija površi.

Parametarske jednačine bilo koje površi zadaju se u obliku:

\[ F(X, Y, Z)=0 \]

Kod sfere imamo slučaj da je rastojanje bilo koje tačke od površi do centra sfere jednako poluprečniku sfere:

\[ F(X, Y, Z)=X^2+Y^2+Z^2-R^2=0 \]

Gde su pravougle koordinate definisane u odnosu na krivolinijiske, kao što je definisano u prethodnom poglavlju:

\[ \begin{aligned} & \mathrm{X}=\mathrm{R} \cos \varphi \cos \lambda \\ & \mathrm{Y}=\mathrm{R} \cos \varphi \sin \lambda \\ & \mathrm{Z}=\mathrm{R} \sin \varphi . \end{aligned} \]

Parametarske linije se mogu definisati u opštem slučaju za sferu. Jednačina paralele definisana je parametarskim koordinatama tako da je parametar \(\varphi\) konstantan, uzmimo da je označen sa \(c1\).

\[ \begin{aligned} & X=X(c 1, \lambda) \\ & Y=Y(c 1, \lambda) \\ & Z=Z(c 1, \lambda) \end{aligned} \]

Jednačina meridijana definisana je parametarskim koordinatama tako da je parametar \(\lambda\) konstantan, uzmimo da je označen sa c2.

\[ \begin{aligned} & X=X(\varphi, c 2) \\ & Y=Y(\varphi, c 2) \\ & Z=Z(\varphi, c 2) \end{aligned} \]

R

Code

# install.packages("plotly")

library(plotly)
# Warning: package 'plotly' was built under R version 4.2.3
# Warning: package 'ggplot2' was built under R version 4.2.3
R= 6377

# koristimo polarne sferne koordinate u je theta, v je fi
theta <- seq(0, 2*pi, length.out = 120)
phi <- seq(0, pi, length.out = 60)

u <- outer(theta, phi, FUN = function(t, p) t)
v <- outer(theta, phi, FUN = function(t, p) p)

# površ sfere
xs <- R*cos(u) * sin(v)
ys <- R*sin(u) * sin(v)
zs <- R*cos(v)

x <- c() ; y <- c(); z <- c()

# kreiranje meridijana gde je jedan parametar konstantan
for (c1 in theta[seq(1, length(theta), by = 10)]) {
  x <- c(x, R*cos(c1) * sin(phi), NA)  # NA za kraj meridijana
  y <- c(y, R*sin(c1) * sin(phi), NA)
  z <- c(z, R*cos(phi), NA)
}

# kreiranje paralela gde je jedan parametar konstantan
for (c2 in phi[seq(1, length(phi), by = 6)]) {
  x <- c(x, R*cos(theta) * sin(c2), NA)  # NA za kraj paralele
  y <- c(y, R*sin(theta) * sin(c2), NA)
  z <- c(z, R*rep(cos(c2), length(theta)), NA)
}

fig <- plot_ly() %>%
  add_surface(x = xs, y = ys, z = zs, colorscale = list(c(0, '#ffffff'), c(1, '#ffffff')),
              showscale = FALSE, opacity = 0.5) %>%
  add_trace(x = x, y = y, z = z, type = 'scatter3d', mode = 'lines',
            line = list(width = 3, color = 'rgb(10,10,10)'))

fig <- fig %>% layout(scene = list(aspectratio = list(x = 1, y = 1, z = 1)))
fig

Slika 3.1: Parametarske linije na sferi.

Python

Code

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plot_sfera(radius=1, meridians=12, parallels=6):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    # polarne sferne koordinate
    u = np.linspace(0, 2 * np.pi, meridians)
    v = np.linspace(0, np.pi, parallels)
    
    x = radius * np.outer(np.cos(u), np.sin(v))
    y = radius * np.outer(np.sin(u), np.sin(v))
    z = radius * np.outer(np.ones(np.size(u)), np.cos(v))

    # sfere
    ax.plot_wireframe(x, y, z, color='b')

    # plot meridijana
    for i in range(meridians):
        ax.plot(x[i, :], y[i, :], z[i, :])

    # plot paralela
    for j in range(parallels):
        ax.plot(x[:, j], y[:, j], z[:, j])

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Parametarske linije na sferi')

    # prikaz
    plt.show()

# plot (# može se napraviti plot i koristeći plotly iz Pythona)
plot_sfera(radius=6377, meridians=12, parallels=6)

Parametarske linije na sferi.

3.2 Parametarske jednačine elipsoida

Parametarske linije na elipsoidu definisane su kao kod sfere, i to su meridijani i paralele na elipsoidu. Jednačina elipsoida ima oblik:

\[ F(X, Y, Z)=X^2+Y^2+Z^2-N^2=0 \]

Gde su:

\[ \begin{aligned} & X=(N+h) \cos \varphi \cos \lambda \\ & Y=(N+h) \cos \varphi \sin \lambda \\ & Z=\left(N\left(1-e^2\right)+h\right) \sin \varphi \end{aligned} \]

Poluprečnik krivine prvog vertikala:

\[ N(\varphi)=\frac{a}{\sqrt{1-e^2 \sin ^2 \varphi}} \]

Prvi numerički ekscentricitet:

\[ e^2=\left(a^2-b^2\right) / a^2 \]

Jednačina paralele definisana je parametarskim koordinatama tako da je parametar \(\varphi\) konstantan, uzmimo da je označen sa \(c1\).

\[ \begin{aligned} & X=X(c 1, \lambda) \\ & Y=Y(c 1, \lambda) \\ & Z=Z(c 1, \lambda) \end{aligned} \]

Jednačina meridijana definisana je parametarskim koordinatama tako da je parametar \(\lambda\) konstantan, uzmimo da je označen sa c2.

\[ \begin{aligned} & X=X(\varphi, c 2) \\ & Y=Y(\varphi, c 2) \\ & Z=Z(\varphi, c 2) \end{aligned} \]

3.3 Gausove fundamentalne veličine

Ako je tačka A definisana koordinatama (X,Y,Z), Slika 3.2. Uočimo tačku C koja je veoma bliska tački A, i definisana je koordinatama X+ΔX, Y+ΔY, Z+ΔZ. Ako smatramo da su A i C vrlo bliske tačke, možemo da zanemarimo zakrivljenost elipsoida i dužinu možemo da sračunamo kao pravu u prostoru.

Slika 3.2: Elementarni (beskonačno mali) elipsoidni trapez..

Onda je dužina AC:

\[ A C^2=(X+\Delta X-X)^2+(Y+\Delta Y-Y)^2+(Z+\Delta Z-Z)^2 \] Odnosno :

\[ A C^2=\Delta X^2+\Delta Y^2+\Delta Z^2 \]

Ako pretpostavimo da su tačke na neizmerno malom rastojanju onda je jasno da veličinu \(ds\) možemo izraziti kao:

\[ d s^2=d X^2+d Y^2+d Z^2 \] Uzmimo da je u pitanju bilo koja površ koja ima parametarske linije \(\varphi\) i \(\lambda\), onda se svaka pravougla koordinata može definisati kao funkcija od parametarskih koordinata:

\[ \begin{aligned} & X=x(\varphi, \lambda) \\ & Y=y(\varphi, \lambda) \\ & Z=z(\varphi, \lambda) \end{aligned} \]

Diferencijali su onda:

\[ \mathrm{dX}=\frac{\partial(x(\varphi, \lambda))}{\partial \varphi} \mathrm{d} \varphi+\frac{\partial(x(\varphi, \lambda))}{\partial \lambda} \mathrm{d} \lambda \]

\[ d Y=\frac{\partial(y(\varphi, \lambda))}{\partial \varphi} d \varphi+\frac{\partial(y(\varphi, \lambda))}{\partial \lambda} d \lambda \]

\[ \mathrm{dZ}=\frac{\partial(z(\varphi, \lambda))}{\partial \varphi} \mathrm{d} \varphi+\frac{\partial(z(\varphi, \lambda))}{\partial \lambda} \mathrm{d} \lambda \]

Ako bi se dX, dY i dZ uvrstile u izraz \(d s^2=d X^2+d Y^2+d Z^2\), dobijamo dugačak izraz koji se skraćeno može zapisati:

\[ d s^2=e d \varphi^2+2 f d \varphi d \lambda+g d \lambda^2 \]

Veličine e, f i g se nazivaju Gausove fundamentalne veličine

\[ \mathrm{e}=\left(\frac{\partial(x(\varphi, \lambda))}{\partial \varphi}\right)^2+\left(\frac{\partial(y(\varphi, \lambda))}{\partial \varphi}\right)^2+\left(\frac{\partial(z(\varphi, \lambda))}{\partial \varphi}\right)^2 \]

\[ f=\frac{\partial(x(\varphi, \lambda))}{\partial \varphi} \frac{\partial(x(\varphi, \lambda))}{\partial \lambda}+\frac{\partial(y(\varphi, \lambda))}{\partial \varphi} \frac{\partial(y(\varphi, \lambda))}{\partial \lambda}+\frac{\partial(z(\varphi, \lambda))}{\partial \varphi} \frac{\partial(z(\varphi, \lambda))}{\partial \lambda} \]

\[ g=\left(\frac{\partial(x(\varphi, \lambda))}{\partial \lambda}\right)^2+\left(\frac{\partial(y(\varphi, \lambda))}{\partial \lambda}\right)^2+\left(\frac{\partial(z(\varphi, \lambda))}{\partial \lambda}\right)^2 \]

Gausove fundamentalne veličine za elipsoid

Za nemerljivo mala rastojanja (rastojanja teže nuli) trougao ADC, Slika 3.3, možemo smatrati ravnim.

Slika 3.3: Beskonačno mali elipsoidni trapez.

Meridijani i paralele se seku pod pravim uglom, tako da je:

\[ \mathrm{ds}^2=(\mathrm{Md} \varphi)^2+(N \cos \varphi \mathrm{d} \lambda)^2 \]

\[ \begin{aligned} & \mathrm{e}=\mathrm{M}^2 \\ & \mathrm{f}=0 \\ & \mathrm{g}=(N \cos \varphi)^2 \end{aligned} \]

3.4 Opšta teorija kartografskog preslikavanja

Za kartografsko preslikavanje između originalne i projekcione površi koriste se neke matematičke funkcije preslikavanja koje preslikavaju tačke, linije i poligone sa elipsoida na ravan. Svakoj tački sa elipsoida odgovara samo jedna tačka u projekciji i obratno. Sa slike, (Slika 3.4), može se uočiti da se tačka A preslikala u tačku A’, tačka B u B’, tačka C u C’ i konačno tačka D u D’. Pored toga, vidimo da u projekciji imamo sliku meridijana i paralela koji su preslikani u opštem slučaju kao krive linije.

Slika 3.4: Ilustracija kartografskog preslikavanja sa elipsoida na ravan karte, sa elipsoidnih krivolinijskih u pravougle koordinate u ravni.

Funkcije kartografskog preslikavanja imaju oblik:

\[E= y = f_{1}(\varphi,\lambda)\]

\[N = x = f_{2}(\varphi,\lambda)\]

Ili za inverzni zadatak:

\[\varphi = g_{1}(E,N)\]

\[\lambda= g_{2}(E,N)\]

Funkcije kartografskog preslikavanja moraju da zadovolje sledeće uslove:

• funkcije preslikavanja su obostrano jednoznačne,

• da su neprekidne sa izvodima 1. i 2. reda,

• Jakobijan različit od nule u svakoj tački.

Formula za Jakobijan u opštem slučaju je:

\[ J_F(M)=\left(\begin{array}{ccc} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \cdots & \frac{\partial f_m}{\partial x_n} \end{array}\right) \]

Linijski element u ravni karte

Sledeća slika je verovatno najznačajnija kod opšte teorije kartografskog preslikavanja zato što mnoge veličine u pogledu preslikavanja sa elipsoida na ravan zaključujemo polazeći od preslikavanja elementarnog tapeza ABCD na ravan karte. Kao što smo ranije pominjali, tačke ovog trapeza su na infinitezimalnom rastojanju, teži beskonačno malom rastojanju, tako da se ovakav trapez može smatrati ravnim. Na slici (Slika 3.5) su prikazane sledeće veličine koje dosad nisu pominjane. Azimut „α“ na elipsoidu utvrđuje položaj pravaca, a meri se počev od meridijana tačke u pravcu kazaljke na časovniku,odnosno od severa na istok, jug i zapad. A je azimut u projekciji, isprekidane linije su tangente na odgovarajuće krive u ravni. Ugao \(\theta\) je ugao između meridijana i paralele u projekciji. Na elipsoidu se meridijani i paralele seku pod pravim uglom, dok se u projekciji seku pod uglom \(\theta\), a samo kod konformnih projekcija pod pravim uglom. Treba napomenuti da je zakrivljenost meridijana i paralela na beskonačno malom trapezu u ravni projekcije zanimarljiva, ali je ovde istaknuta zbog razumevanja procesa kartografskog preslikavanja.

Slika 3.5: Preslikavanje infinitezimalnog trapeza sa elipsoida na ravan karte.

Vratimo se na opšte jednačine kartografskog preslikavanja:

\[ \begin{aligned} & y=\text { Eeasting }=f_1(\varphi, \lambda)=y(\varphi, \lambda) \\ & x=\text { Northing }=f_2(\varphi, \lambda)=x(\varphi, \lambda) \end{aligned} \]

Dužina dS u ravni može se sračunati koristeći Gausove fundamentalne veličine, već smo napomenuli i ranije da one mogu da važe za bilo koju površ, pa i ravan u kojoj su pravougle koordinate funkcije od krivolinijskih.

\[ d\mathrm{S}^2=\mathrm{Ed} \varphi^2+2 \mathrm{Fd} \varphi \mathrm{d} \lambda+\mathrm{Gd} \lambda^2 \]

Gde su E, F i G Gausove fundamentalne veličine za ravan, a pravougle koordinate su funkcije od krivolinijskih elipsoidnih koordinata.

\[ \mathrm{E}=\left(\frac{\partial(x(\varphi, \lambda))}{\partial \varphi}\right)^2+\left(\frac{\partial(y(\varphi, \lambda))}{\partial \varphi}\right)^2 \]

\[ \mathrm{F}=\frac{\partial(x(\varphi, \lambda))}{\partial \varphi} \frac{\partial(x(\varphi, \lambda))}{\partial \lambda}+\frac{\partial(y(\varphi, \lambda))}{\partial \varphi} \frac{\partial(y(\varphi, \lambda))}{\partial \lambda} \]

\[ \mathrm{G}=\left(\frac{\partial(x(\varphi, \lambda))}{\partial \lambda}\right)^2+\left(\frac{\partial(y(\varphi, \lambda))}{\partial \lambda}\right)^2 \]

Znači da sada za bilo koju kartografsku projekciju koja je zadata u obliku:

\[ y=\text { Easting }=f_1(\varphi, \lambda)=y(\varphi, \lambda)\] ,
\[ x=\text { Northing }=f_2(\varphi, \lambda)=x(\varphi, \lambda) \]

možemo sračunati linijski element dS. Još jednom treba napomenuti da je kod konformnih projekcija, gde se meridijani i paralele seku pod pravim uglom, veličina F jednaka nuli.

Linearna razmera

Linearnim razmera je granična vrednost (limes) odnosa beskonačno male duži u projekciji i odgovarajuće duži na površi Zemljinog elipsoida ili lopte. Stoga je, Slika 3.5, formula za linearnu razmeru:

\[ \mathrm{c}=\lim _{A C \rightarrow 0}\left(\frac{A^{\prime} C^{\prime}}{A C}\right)=\frac{d S}{d s} \] Ako kvadriramo levu i desnu stranu, koristeći Gausove fundamentalne veličine, imamo:

\[ c^2=\frac{E d \varphi^2+2 \mathrm{Fd} \varphi \mathrm{d} \lambda+\mathrm{Gd} \lambda^2}{(\mathrm{Md} \varphi)^2+(N \cos \varphi \mathrm{d} \lambda)^2} \]

Uočimo sa slike (Slika 3.6) da azimut pravca \(ds\) možemo izraziti koristeći dužinu luka meridijana i paralele kod infinitezimalnog trapeza.

Slika 3.6: Infinitezimalni trapez na elipsoidu.

Tako da je:

\[ \tan \alpha=\frac{N \cos \varphi \mathrm{d} \lambda}{\operatorname{Md} \varphi} \]

Odatle sledi da je:

\[ \frac{\mathrm{d} \lambda}{\mathrm{d} \varphi}=\frac{\mathrm{M} \tan \alpha}{N \cos \varphi} \]

Dalje, ako jednačinu za kvadrat linearne razmere pomnožimo i podelimo sa \(1/\mathrm{d}\varphi^2\) tada imamo:

\[ \mathrm{c}^2=\frac{\mathrm{Ed} \varphi^2+2 \mathrm{Fd} \varphi \mathrm{d} \lambda+\mathrm{Gd} \lambda^2}{(\mathrm{Md} \varphi)^2+(N \cos \varphi \mathrm{d} \lambda)^2} \frac{\frac{1}{\mathrm{~d} \varphi^2}}{\frac{1}{\mathrm{~d} \varphi^2}}=\frac{\mathrm{E}+2 \mathrm{~F} \frac{\mathrm{d} \lambda}{\mathrm{d} \varphi}+\mathrm{G}\left(\frac{\mathrm{d} \lambda}{\mathrm{d} \varphi}\right)^2}{\mathrm{M}^2+(N \cos \varphi)^2\left(\frac{\mathrm{d} \lambda}{\mathrm{d} \varphi}\right)^2} \]

Zamenimo \(\frac{\mathrm{d} \lambda}{\mathrm{d} \varphi}\) i dobijamo:

\[ \mathrm{c}^2=\frac{\mathrm{E}+2 \mathrm{~F} \frac{\mathrm{M} \tan \alpha}{N \cos \varphi}+\mathrm{G}\left(\frac{\mathrm{M} \tan \alpha}{N \cos \varphi}\right)^2}{\mathrm{M}^2+(N \cos \varphi)^2\left(\frac{\mathrm{M} \tan \alpha}{N \cos \varphi}\right)^2} \]

ili:

\[ \mathrm{c}^2=\frac{\mathrm{E}+2 \mathrm{~F} \frac{\mathrm{M} \tan \alpha}{N \operatorname{cos} \varphi}+\mathrm{G}\left(\frac{\mathrm{M} \tan \alpha}{N \cos \varphi}\right)^2}{\mathrm{M}^2\left(1+\tan ^2 \alpha\right)} \]

Već odavde je jasno da je linearna razmera funkcija koja zavisi od pozicije tačke i azimuta pravca

\[ c=c(\varphi, \lambda, \alpha) \]

Dakle, u opštem slučaju na svakoj poziciji na karti linearna razmera ima drugu vrednost. Na samoj poziciji u različitim pravcima ima različite vrednosti.

Primenom trigonometrijske formule:

\[ \sin (2 \alpha)=2 \sin \alpha \cos \alpha= (\sin \alpha +\cos \alpha)^2 - 1 =\frac{2 \tan \alpha}{1+\tan ^2 \alpha} \]

i sređivanjem izraza došlo bi se do konačne formule za linearnu razmeru u opštem slučaju.

\[ \mathrm{c}^2=\frac{E}{M^2} \cos ^2 \alpha+\frac{F}{M N \cos \varphi} \sin 2 \alpha+\frac{G}{(N \cos \varphi)^2} \sin ^2 \alpha \] Odstupanje razmere od 1 naziva se linearnom deformacijom ili deformacijom dužina (\(dc\)):

\[ dc=c-1 \]

Linearna razmera u pravcu meridijana i paralele

Uočimo sa slike (Slika 3.5) da je linearna razmera u pravcu meridijana (m):

\[ m=\frac{\sqrt{E}}{M} \]

Ovo smo mogli dobiti i računajući preko opšte formule, azimut se meri od meridijana, jasno je da je od tačke A u pravcu meridijana azimut 0, Slika 3.5.

\[ \mathrm{c}_{\mathrm{m}}^2=m^2=\frac{E}{M^2} \cos ^2 0+\frac{F}{M N \cos \varphi} \sin 0+\frac{G}{(N \cos \varphi)^2} \sin ^2 0 \]

Slično sa slike (Slika 3.5) jasno je da je linearna razmera u pravcu paralele:

\[ n=\frac{\sqrt{G}}{N \cos \varphi} \]

Ili, ako bismo zamenili u formuli za azimut 90 stepeni imali bismo izvođenje za linearnu razmeru u pravcu paralele:

\[ \mathrm{c}_{\mathrm{n}}{ }^2=n^2=\frac{E}{M^2} \cos ^2 90+\frac{F}{M N \cos \varphi} \sin 180+\frac{G}{(N \cos \varphi)^2} \sin ^2 90 \]

Linearna razmera kao funkcija m i n

Ako uočimo sliku trapeza na elipsoidu koji teži da bude beskonačno mali, njegova slika na karti je beskonačno mali paralelogram, gde za tako mali element slike meridijana i paralele možemo smatrati pravim linijama u okviru paralelograma, Slika 3.7.

Slika 3.7: Infinitezimalni trapez na elipsoidu i njegova slika u ravni za izvođenje formule za linearnu razmeru u funkciji m i n.

Tada linijski element možemo sračunati koristeći kosinusnu teoremu i onda imamo:

\[ \mathrm{d} \mathrm{S}^2=\mathrm{Ed} \varphi^2-2 \sqrt{E G} \cos (180-\theta)+\mathrm{G} \mathrm{d} \lambda^2 \]

dalje je:

\[ \mathrm{dS}{ }^2=\mathrm{Ed} \varphi^2+2 \sqrt{E G} \cos (\theta)+\mathrm{G} \mathrm{d} \lambda^2 \]

Ako se vratimo na isti element izražen koristeći Gausove fundamentalne veličine:

\[ \mathrm{d} \mathrm{S}^2=\mathrm{Ed} \varphi^2+2 \mathrm{Fd} \varphi \mathrm{d} \lambda+\mathrm{Gd} \lambda^2 \]

Možemo videti da su kvadratni članovi isti, tako da mora biti isti i mešoviti član, pa dalje imamo da je kosinus ugla između meridijana i paralele u projekciji:

\[ \cos (\theta)=\frac{\mathrm{F}}{\sqrt{E G}} \]

Odavde dalje možemo sračunati i F kao:

\[ F=\sqrt{E G} \cos (\theta) \]

Zamenimo F u opštoj formuli za linearnu razmeru \(\mathrm{c}^2=\frac{E}{M^2} \cos ^2 \alpha+\frac{F}{M N \cos \varphi} \sin 2 \alpha+\frac{G}{(N \cos \varphi)^2} \sin ^2 \alpha\) i dobijamo:

\[ \mathrm{c}^2=\frac{E}{M^2} \cos ^2 \alpha+\frac{\sqrt{E G} \cos (\theta)}{M N \cos \varphi} \sin 2 \alpha+\frac{G}{(N \cos \varphi)^2} \sin ^2 \alpha \]

Ako se vratimo na formulu za linearnu razmeru u pravcu meridijana:

\[ m=\frac{\sqrt{E}}{M} \]

i paralele:

\[ n=\frac{\sqrt{G}}{N \cos \varphi} \]

Tada dobijamo konačnu formulu za linearnu razmeru u finkciji linearne razmere u pravcu meridijana i paralele:

\[ \mathrm{c}^2=m^2 \cos ^2 \alpha+m n \cos (\theta) \sin 2 \alpha+n^2 \sin ^2 \alpha \]

Glavni pravci

Pravci u kojima linearna razmera dostiže ekstremnu vrednost nazivaju se glavni pravci

\[ \mathrm{c}^2=\frac{E}{M^2} \cos ^2 \alpha+\frac{F}{M N \cos \varphi} \sin 2 \alpha+\frac{G}{(N \cos \varphi)^2} \sin ^2 \alpha \]

Ako uvedemo zamenu (R,P,Q) za elemente koji ne zavise od azimuta, formula je:

\[ c^2=R \cos ^2 \alpha+P \sin 2 \alpha+Q \sin ^2 \alpha \]

\[ \frac{\partial c}{\partial \alpha}=-2 \mathrm{R} \cos \alpha \sin \alpha+2 \mathrm{P} \cos 2 \alpha+2 \mathrm{Q} \sin \alpha \cos \alpha \]

\[ \frac{\partial c}{\partial \alpha}=2 \mathrm{P} \cos 2 \alpha+2(\mathrm{Q}-\mathrm{R}) \sin \alpha \cos \alpha \]

\[ \frac{\partial c}{\partial \alpha}=2 \mathrm{P} \cos 2 \alpha+(\mathrm{Q}-\mathrm{R}) \sin 2 \alpha \]

Ekstremne vrednosti linearne razmere su u pravcima u kojim je prvi izvod 0.

\[ \frac{\partial c}{\partial \alpha}=2 \mathrm{P} \cos 2 \alpha_0+(\mathrm{Q}-\mathrm{R}) \sin 2 \alpha_0=0 \]

Ako izraz podelimo sa \(\cos 2 \alpha_0\)

\[ \begin{aligned} & 2 \mathrm{P}+(\mathrm{Q}-\mathrm{R}) \tan 2 \alpha_0=0 \\ & \tan 2 \alpha_0=\frac{2 P}{R-Q} \end{aligned} \]

Funkcija tangensa

Slika 3.8: Tangensna funkcija, (https://sh.wikipedia.org/wiki/Tangens#/media/Datoteka:Tan.svg).

Pošto je period tangensa jednak \(180^o\) , jednačina ima 2 rešenja: jedno će biti \(\alpha_0\), a drugo \(\alpha_0+90\). Prema tome, glavni pravci se na originalnoj površi seku pod uglom od \(90^o\).

Rešenja su oblika: \[ \alpha_0 \text { i } \alpha_0+90^o \]

Glavni pravci su upravni jedan na drugi, seku se pod pravim uglom, i u projekciji (ravni karte) i na površi (elipsoidu ili lopti) su ortogonalni (upravni). Zaključak je da: uvek postoje samo dva upravna pravca na elipsoidu koji će se preslikati na karti kao ortogonalni (upravni) i imaće ekstremnu vrednost linearne razmere. Pravac u kom je linearna razmera najveća naziva se prvi glavni pravac i obično se obeležava slovom a. Pravac u kom je linearna razmera najmanja naziva se drugi glavni pravac i obično se obeležava slovom b.

Razmer površina

Odnos površina odgovarajućih elementarnih figura na slici i originalnoj površi naziva se razmer površina, Slika 3.7. Površina beskonačno malog trapeza na elipsoidu je:

\[ \mathrm{dp}=\mathrm{MN} \cos \varphi \mathrm{d} \lambda \mathrm{d} \varphi \]

U ravni karte ako posmatramo slike: Slika 3.7, Slika 3.9,

Slika 3.9: Površina elemantarnog paralelograma u ravni karte.

možemo zaključiti da je površina odgovarajućeg paralelograma:

\[ \mathrm{dP}=\sqrt{G} \mathrm{~d} \lambda \sqrt{E} \mathrm{~d} \varphi \sin (\theta) \]

U delu kod računanja linerane razmere u funkciji m i n smo videli da je:

\[ \cos (\theta)=\frac{\mathrm{F}}{\sqrt{E G}} \]

Koristeći trigonometrijske osobine sinusa i kosinusa imamo:

\[ \begin{aligned} & \sin ^2(\theta)=1-\cos ^2(\theta) \\ & \sin ^2(\theta)=1-\frac{F^2}{E G}=\frac{E G-F^2}{E G} \\ & \sin (\theta)=\sqrt{\frac{E G-F^2}{E G}}=\frac{H}{\sqrt{E G}} \\ & \end{aligned} \]

Tako da imamo dalje, ako sa H ozačimo \(H= \sqrt{E G-F^2}\), da je:

\[ \mathrm{dP}=\sqrt{G} \mathrm{~d} \lambda \sqrt{E} \mathrm{~d} \varphi \sin (\theta)=\sqrt{E G} \frac{H}{\sqrt{E G}} \mathrm{~d} \lambda \mathrm{d} \varphi \]

I konačno imamo:

\[ \mathrm{dP}=H \mathrm{~d} \lambda \mathrm{d} \varphi \]

Formula za razmer površina je:

\[ p=\frac{\mathrm{dP}}{\mathrm{dp}}=\frac{\mathrm{Hd} \lambda \mathrm{d} \varphi}{\mathrm{MN} \cos \varphi \mathrm{d} \lambda \mathrm{d} \varphi}=\frac{\mathrm{H}}{\mathrm{MN} \cos \varphi} \]

Pošto je \(\sin (\theta)=\frac{H}{\sqrt{E G}}\), onda je i \(\mathrm{H}=\sqrt{E G} \sin (\theta)\) tada formula ima oblik:

\[ p=\frac{\sqrt{E G} \sin (\theta) \mathrm{d} \lambda \mathrm{d} \varphi}{\mathrm{MN} \cos \varphi \mathrm{d} \lambda \mathrm{d} \varphi}=\mathrm{mn} \sin (\theta) \]

Ranije smo videli da je \(m=\frac{\sqrt{E}}{M}\) i \(n=\frac{\sqrt{G}}{N \cos \varphi}\).

Deformacija površina je definisana izrazom koji pokazuje koliko se one deformišu u odnosu na preslikavanje koje nema deformacija površina, tj. gde je p=1:

\[ d p=p-1 \]

Deformacije uglova

U procesu kartografskog preslikavanja ne deformišu se samo dužine i površine nego i uglovi. Da bi izvođenje bilo jednostavnije, pretpostavićemo da se koordinatne ose poklapaju sa glavnim pravcima (izvođenje se može primeniti i u opštem slučaju). Uočimo beskonačno mali krug na elipsoidu, Slika 3.10, postavimo pravougli koordinatni sistem X,Y čija je slika pravougli koordinatni sistem u ravni koji se poklapa sa glavnim pravcima a i b, a koordinatne ose su označene kao \(x,y\). Imamo da je slika beskonačno malog kruga elipsa (elipsa deformacija ili Tisoova idikatrisa). Videli smo da je u opštem slučaja linearna razmera različita u različitim pravcima, imaju dva upravna glavna pravca u kojim ona dostiže ekstremne vrednosti a i b, razmera u pravcu paralele je n i u pravcu meridijana je m. Azimut glavnog pravca na elipsoidu je označen sa \(\alpha\), dok je ugao \(\beta\) direkcioni ugao tačke T sa elementarnog kruga u odnosu na koordinatni sistem postavljen na elipsoidu (\(X,Y\)). Tačka T se preslikava u tačku T’ koja ima direkcioni ugao B, a azimut glavnog pravca je označen sa A u ravni projekcije.

Slika 3.10: Deformacija uglova.

Deformacija ugla predstavlja razliku ugla na elipsoidu i ugla u ravni (projekciji), obeležava se sa \(du\) ili \(\omega/2\)

\[ \omega / 2=\beta-B \]

Sa slike možemo odrediti glavne pravce kao:

\[ a=\frac{x}{X}, b=\frac{y}{Y} \]

Zato što su koordinate tačke T u pravcu \(x\) i \(y\) infinitezimalne dužine u pravcu glavnih pravaca, koje su slika odgovarajućih \(X\) i \(Y\) dužina sa elipsoida. Dalje imamo:

\[ \frac{y}{x}=\frac{b Y}{a X} \]

Pošto je \(\frac{y}{x}\) u stvari nagib pravca tačke T’, a \(\frac{Y}{X}\) nagib pravca tačke T, imamo:

\[ \tan B =\frac{b}{a} \tan \beta \]

Uočimo izraz

\[ \frac{\sin (\beta-B)}{\sin (\beta+B)}=\frac{\sin \beta \cos B-\cos \beta \sin B}{\sin \beta \cos B+\cos \beta \sin B} \]

Podelimo brojilac i imenilac sa desne strane izraza sa \(cos \beta cos B\):

\[ \frac{\sin (\beta-B)}{\sin (\beta+B)}=\frac{\operatorname{tan} \beta-\frac{b}{a} \operatorname{tan} \beta}{\operatorname{tan} \beta+\frac{b}{a} \operatorname{tan} \beta}=\frac{a-b}{a+b} \]

Odakle konačno dobijamo:

\[ \sin (\beta-B)=\frac{a-b}{a+b} \sin (B+\beta) \]

Nas zanima maksimalna deformacija uglova u okolini posmatrane tačke. Pošto a i b imaju konstantne vrednosti u zadatoj tački karte, očigledno je da će razlika \((\beta-B)=du=\omega/2\) dostići maksimalan iznos kada je \(sin (\beta+B)=1\), odnosno kada je zbir uglova \((\beta+B)=90^o\).

\[ \sin (\beta-B)=\frac{a-b}{a+b} \]

\[ \sin \frac{\omega}{2}=\frac{a-b}{a+b} \]

U litaraturi, npr. Jovanović (1984), mogu se naći i sledeći izrazi za maksimalnu deformaciju uglova:

\[ \cos \frac{\omega}{2}=\sqrt{1-\sin ^2 \frac{\omega}{2}}=\frac{2 \sqrt{a b}}{a+b} \] \[\operatorname{tan} \frac{\omega}{2}=\frac{\sin \frac{\omega}{2}}{\cos \frac{\omega}{2}}=\frac{a-b}{2 \sqrt{a b}}\]

\[ \operatorname{tan} \frac{\omega}{4}=\sqrt{\frac{1-\cos \frac{\omega}{2}}{1+\cos \frac{\omega}{2}}}=\frac{\sqrt{a}-\sqrt{b}}{\sqrt{a}-\sqrt{b}} \]

\[ \operatorname{tan}\left(45^0+\frac{\omega}{4}\right)=\frac{1+\operatorname{tan} \frac{\omega}{4}}{1-\operatorname{tan} \frac{\omega}{4}}=\sqrt{\frac{a}{b}} \]

\[ \operatorname{tan}\left(45^{\circ}-\frac{\omega}{4}\right)=\frac{1-\operatorname{tan} \frac{\omega}{4}}{1+\operatorname{tan} \frac{\omega}{4}}=\sqrt{\frac{a}{b}} \]

\[ \sin \frac{\omega}{2}=\frac{a-b}{a+b}=\sqrt{\frac{m^2+n^2-2 m n \sin \theta}{m^2+n^2+2 m n \sin \theta}}=\sqrt{\frac{m^2+n^2-2 m n \cos \varepsilon}{m^2+n^2+2 m n \cos \varepsilon}} \]

\[ \operatorname{tan} \frac{\omega}{2}=\frac{a-b}{2 \sqrt{a b}}=\frac{a-b}{2 \sqrt{p}} \]

\[ \operatorname{tan} \frac{\omega}{2}=\frac{1}{2} \sqrt{\frac{m^2+n^2-2 m n \sin \theta}{m n \sin \theta}}=\frac{1}{2} \sqrt{\frac{m^2+n^2-2 m n \cos \varepsilon}{m n \cos \varepsilon}} \]

Gde je \(\varepsilon = \theta -90\), odstupanje ugla između meridijana i paralele u projekciji od pravog ugla.

3.5 Elipsa deformacija (Tisoova indikatrisa)

Francuski geograf i matematičar Tiso, razmatrajući deformacije koje nastaju pri preslikavanju dvostruko zakrivljene površi Zemlje na ravan, došao je do veoma značajnih zaključaka za matematičku kartografiju.

Teorija deformacija po Tisou, u osnovi, objašnjava da se elementarni krug sa površi elipsoida, načelno uzimajući, preslikava na ravan kao elipsa, osim u specijalnom slučaju ako su obe površi (ili njihovi delovi) međusobno paralelne, u kom slučaju će krug ostati krug željenog razmera.

Razmotrimo slučaj da se postavi dvodimenzionalni koordinatni sistem na elipsoidu tako da se poklapa sa pravcima meridijana i paralele, i da se iz ishodišta kreira infinitezimalni krug, Slika 3.11.

Slika 3.11: Infinitezimalni krug na elipsoidu i odgovarajuće preslikavanje u ravni pravca OT.

Možemo uočiti element \(ds\):

\[ \mathrm{d} X^2+\mathrm{d} Y^2=\mathrm{d} s^2 \]

Linearna razmera u pravcu meridijana:

\[ \frac{O^{\prime} \cdot P^{\prime}}{O P}=\frac{\mathrm{d} x}{\mathrm{~d} X}=m \]

Linearna razmera u pravcu paralele:

\[ \frac{O^{\prime} \cdot Q^{\prime}}{O Q}=\frac{\mathrm{d} y}{\mathrm{~d} Y}=n \]

Ako linijski element na elipsoidu izrazimo koristeći koordinate iz ravni projekcije imamo:

\[ \frac{\mathrm{d} x^2}{m^2}+\frac{\mathrm{d} y^2}{n^2}=\mathrm{d} s^2 \]

Potom izraz možemo podeliti sa \(\mathrm{d} s^2\) i dobijammo:

\[ \frac{\mathrm{d} x^2}{m^2 \mathrm{~d} s^2}+\frac{\mathrm{d} y^2}{n^2 \mathrm{~d} s^2}=1 \]

Uočavamo da navedena formula predstavlja formulu elipse, pri čemu su \(m ds\) i \(n ds\) tzv. spregnuti (konjugovani) dijametri elipse. Svaki par međusobno upravnih dijametara kruga preslikava kao par konjugovanih dijametara elipse. U opštem slučaju vidimo da se elemantarni krug se preslikava u elipsu, Tisoovu indikatrisu. Tisoova indikatrisa nam pokazuje karakteristike deformacija u okolini posmatrane tačke, sa nje uočavamo veličine deformacija i to: deformacije uglova (oblika), dužina i površina.

Uzmimo da je ds=1, neka jedinična vrednost, tada imamo:

\[ \frac{\mathrm{d} x^2}{m^2}+\frac{\mathrm{d} y^2}{n^2}=1 \]

Ako su koordinate (postavljene kroz centar elementarnog kruga, Slika 3.12) na elipsoidu ujedno i glavni pravci, drugim rečima, ako mreža meridijana i paralela ostaje ortogonalna i u projekciji, onda će, uz uslov da se prvi glavni pravac poklapa sa meridijanom, biti: m=a i n=b, i tada sledi da je:

\[ \frac{\mathrm{d} x^2}{a^2}+\frac{\mathrm{d} y^2}{b^2}=1 \]

Slika 3.12: Elipsa deformacija kada se glavni pravci poklapaju sa koordinatnim sistemom definisanim u centru infinitezimalnog kruga.

Apolonijeve teoreme nam mogu pomoći pri određivanju vrednosti prvog i drugog glavnog pravca (\(a\) i \(b\)) i konstruisanju elipse deformacija. Prema prvoj teoremi, suma kvadrata spregnutih dijametara elipse jednaka je sumi kvadrata njenih poluosa:

\[ a^2+b^2=m^2+n^2 \]

Druga Apolonijeva teorema glasi: površina tangentnog paralelograma elipse nad njenim dijametrima jednaka je konstanti, odnosno površini tangentnog paralelograma nad poluosama elipse, tj.:

\[ a \cdot b=m \cdot n \cdot \sin \theta \]

Na osnovu Apolonijevih teorema imamo:

\[ \begin{aligned} & a^2+b^2=m^2+n^2 \\ & a \cdot b=m \cdot n \cdot \sin \theta \end{aligned} \]

Ako drugi izraz pomnožimo sa minus dva i saberemo sa prvim dobićemo:

\[ \begin{aligned} &(a+b)^2=m^2+n^2+2 m n \sin \theta\\ \end{aligned} \]

Ako drugi izraz pomonožimo sa minus dva i saberemo sa prvim dobićemo:

\[ \begin{aligned} &(a-b)^2=m^2+n^2-2 m n \sin \theta \end{aligned} \]

Jasno je da onda imamo:

\[ \begin{aligned} &A=(a+b)=\sqrt{m^2+n^2+2 m n \sin \theta}\\ &B=(a-b)=\sqrt{m^2+n^2-2 m n \sin \theta} \end{aligned} \]

Dakle veličine A i B možemo sračunati ako znamo linearnu razmeru u pravcu paralele, u pravcu meridijana i ugao između meridijana i paralele, sve ove veličine, tj. njihove formule, ranije smo izveli. Dalje, ako znamo vrednost za A i B, možemo sračunati vrednosti za prvi i drugi glavni pravac:

\[ \begin{aligned} &a=\frac{A+B}{2}\\ &b=\frac{A-B}{2} \end{aligned} \]

Ako je mreža meridijana i paralela u projekciji ortogonalna, onda je \(\theta= 90^o\), tada je:

\[ \begin{aligned} &a+\mathrm{b}=\mathrm{m}+\mathrm{n}\\ &a-\mathrm{b}=\mathrm{m}-\mathrm{n} \end{aligned} \]

Imamo da je \(a=m\), \(b= n\), ili obrnuto.

Kod konformnih projekcija mreža meridijana i paralela je ortogonalna, jer je ortogonalna i na origanalnoj površi kod sfere ili elipsoida (nema deformacija uglova). Kod konformnih projekcija elipsa deformacija je krug jer nema deformacija oblika. Uslov konformnosti možemo formulisati tako da je u okolini posmatrane tačke linearna razmera u svim pravcima ista.

Elipsa deformacija konstrukcija

Slika 3.13 prikazuje postupak konstrukcije elipse deformacija, slika je podeljena na 5 koraka da bi se lakše pratio tok konstrukcije. Ako znamo linearnu razmeru u pravcu meridijana (\(m\)) i paralele (\(n\)), kao i ugao između meridijana i paralele (\(\theta\)), možemo konstruisati elipsu deformacija.

Slika 3.13: Konstrukcija elipse deformacija.

Dokaz za postupak za grafičku konstrukciju elipse deformacija započinjemo kosinusnom teoremom za trougao CON, Slika 3.14.

Slika 3.14: Konstrukcija elipse deformacija, trougao CON.

Kosinusna teorema za trougao CON

\[ C N^2=m^2+n^2-2 m n \cos (90-\theta) \]

\[ C N^2=m^2+n^2-2 m n \sin (\theta) \]

Kosinusna teorema za trougao DON

\[ D N^2=m^2+n^2-2 m n \cos (90+\theta) \]

\[ D N^2=m^2+n^2+2 m n \sin (\theta) \]

Ako se vratimo na Apolonijevu teoremu i uporedimo sa dobijenim

\[ A=(a+b)=\mathrm{DN}=\sqrt{m^2+n^2+2 m n \sin \theta} \]

\[ B=(a-b)=\mathrm{CN}=\sqrt{m^2+n^2-2 m n \sin \theta} \]

Uočimo da je trougao DCN sličan trouglu DOE, paralelni kraci i pripadajući uglovi, očigledno da je OD = DC/2, tako da važi:

\[ \begin{aligned} & \mathrm{OE}=\mathrm{CN} / 2=(a-b) / 2 \\ & \mathrm{DE}=\mathrm{DN} / 2=(a+b) / 2 \end{aligned} \]

Treba zapaziti da je \(EK = EL = EO\) jer je to poluprečnik malog kruga sa centrom u E, poluprečnik kruga je \(EO\). Tačke \(L\) i \(K\) su dobijene presekom tog kruga sa pravom \(DN\). Pravac velikih poluosa je definisan pravcima \(OL\) i \(OK\). Grafički se mogu preuzeti i vrednosti za veliku i malu poluosu elipse:

\[ \mathrm{DK}=\mathrm{DE}+\mathrm{EL}(\text { isto je kao OE) }=(\mathrm{a}+\mathrm{b}) / 2+(\mathrm{a}-\mathrm{b}) / 2=\mathrm{a} \]

\[ D L=D E-E L(\text { isto je kao OE })=(a+b) / 2-(a-b) / 2=b \]

Pošto elipse deformacija pokazuju prostiranje deformacija u bilo kojoj tački od interesa vrlo se često koriste da sagledamo prostiranje deformacija kod neke projekcije. Besplatan softver flexprojector može se koristiti za pregled rasporeda deformacija kod raznih projekcija koje se koriste za globalno kartiranje planete Zemlje, Slika 3.15.

Slika 3.15: Sinusoidna ekvivalenta projekcija sa Tisoovim indikatrisama.

Elipsa deformacija alternativni način

Kao što smo videli u predhodnom poglavlju R softversko okruženje u okviru paketa s2 i sf ima implementirane algoritme za računanja na elipsoidu ili sferi, isto tako u okviru implementacije PROJ biblioteke moguće je izvršavati direktni i inverzni zadatak, kao i koordinatne transformacije.

Dakle, ako kreiramo konačne krugove na elipsoidu (Slika 3.16) i onda ih transformišemo u neku projekciju, dobićemo elipse deformacija. Ovde treba napomenuti da strogo gledano postupak nije ispravan jer ne koristimo krugove čiji poluprečnik teži nuli, ali suštinski dobijamo prostiranje deformacija na skoro isti način kao da kostruišemo elipse deformacija na gore opisan način (Gimond 2022).

library(sf)

# kreiramo tačke u kojima cemo prikazati elipse deformacija
lonlat = rbind(c(19,42), c(19,44), c(19,46),
               c(20,42), c(20,44), c(20,46),
                c(21,42), c(21,44), c(21,46),
               c(22,42), c(22,44), c(22,46) )
tissot.pt <- st_sfc( st_multipoint(lonlat),  
                     crs = 4326) # WGS84 koordinatni ref. sis.

tissot.pt <- st_cast(tissot.pt, "POINT") 
# konvertujemo u single point geometriju
# plot(tissot.pt)

# koristi se s2 paket za elipsodna računanja
# Pravimo tacke po krugu udaljene 30 km od centra elipse
tissot.sf = st_buffer(tissot.pt, dist = 30000, max_cells = 1000)
# max_cells
# The maximum number of cells to approximate a buffer.
# nrow(st_coordinates(tissot.sf[1]))
# [1] 646
# Ako koristimo opciju max_cell = 1000 imamo 646 prelomnih tačaka poligona

# konvertujemo kolekciju poligona u sf objekat
tissot.sf <-  st_sf(tissot.sf)

# Sračunajmo površinu svakog kruga na elipsoidu
tissot.sf$a <- st_area(tissot.sf)

plotKML::plotKML(tissot.sf)
# [1] 0

Slika 3.16: Krugovi poluprečnika 30 km na WGS84.

Ako transformišemo krugove u Merkatorovu projekciju, dobićemo elipse deformacija u Merkatorovoj projekciji, za svaku elipsu izrazićemo i razmer površina. Uz pomenuti prikaz ćemo kartirati i granice opština u Srbiji.

tissot.merc <- st_transform(tissot.sf, "+proj=merc +ellps=WGS84")
# površine u projekciji
tissot.merc$a_proj <- st_area(tissot.merc)
# odnos površina u projekciji i na elipsoidu
tissot.merc$p <- tissot.merc$a_proj/tissot.merc$a

# granice opština u Srb
if (file.exists('data/pop.gml') ){
    pop = st_read('data/pop.gml')   
  }else{
    url=paste('http://osgl.grf.bg.ac.rs/geoserver/osgl_4/wfs',
          '?service=WFS&version=1.0.0&',
          'request=GetFeature&typeName=osgl_4:POP_age',sep='')

    download.file(url, 'pop.gml')
    pop = st_read('pop.gml')
}
# Reading layer `POP_age' from data source `C:\mk_knjiga\data\pop.gml' using driver `GML'
# Simple feature collection with 197 features and 14 fields
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 18.8 ymin: 41.9 xmax: 23 ymax: 46.2
# CRS:           NA

# podaci su u WGS84
st_crs(pop) =  4326
pop = st_transform(pop, "+proj=merc +ellps=WGS84")


ll = st_centroid(tissot.merc)
ll = st_coordinates(ll)

# Karta sa elipsama deformacija i razmerom površina

plot(st_geometry(pop),  graticule = TRUE, border = 'grey',  axes = TRUE)
plot(st_geometry(tissot.merc),  border = 'red', add = TRUE)
text(ll[, 1], ll[, 2], paste("p=",round(tissot.merc$p,3)) ,  cex = 0.8)

Slika 3.17: Elipse deformacija Merkatorova projekcija.

Kao što vidimo sa slike (Slika 3.17) u Merkatorovoj projekciji površine na severu Srbije deformišu se duplo. U UTM projekciji koja je takođe konformna, površine se znatno manje deformišu, Slika 3.18, to je i razlog zašto se koristi kao državna projekcija u skoro svim zemljama Evrope. Još jednom treba napomenuti da je razmer površina sračunat približno koristeći odnose površina konačnih figura. Razmere površina su sračunate po formulama koje ćemo razraditi u poglavlju koje se odnosi na UTM projekciju.

tissot.utm <- st_transform(tissot.sf, "+proj=utm +zone=34 +ellps=GRS80")
# površine u projekciji
tissot.utm$a_proj <- st_area(tissot.utm)
# računanje linearne razmere u UTM
lr <- function(lonlatmatrica){
    l= (lonlatmatrica[,1] -21)*pi/180
    fi= (lonlatmatrica[,2])*pi/180
    #  GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
    a=6378137.0 
    rf=298.257222101
    f = 1/rf
    e = sqrt(2*f - f^2)
    e1 = sqrt( (2*f - f^2)*(1 - f)^-2 )  
    t=tan(fi)
    eta=e1*cos(fi)
    C1=1/2*(1+eta^2)*(cos(fi))^2
    C2=1/24*(5-4*t^2)*(cos(fi))^4
    # linerarna razmera (u R-u se uglavnom ne koristi c kao promenljiva)
    cc=0.9996*(1+C1*l^2+C2*l^4) 
    return(cc)
}

# linerarna deformacija u dm/km
tissot.utm$p <-  (lr(lonlat)) ^2
# tissot.utm$p <- tissot.utm$a_proj/tissot.utm$a

pop = st_transform(pop, "+proj=utm +zone=34 +ellps=GRS80")

ll = st_centroid(tissot.utm)
ll = st_coordinates(ll)

# Karta sa elipsama deformacija i razmerom površina

plot(st_geometry(pop),  graticule = TRUE, border = 'grey',  axes = TRUE)
plot(st_geometry(tissot.utm),  border = 'red', add = TRUE)
text(ll[, 1], ll[, 2], paste("p=",round(tissot.utm$p,4)) ,  cex = 0.8)

Slika 3.18: Elipse deformacija UTM projekcija.

Kod Sinusoidne projekcije površine se ne deformišu jer je ona ekvivalentna, ali se deformišu oblici, Slika 9.4. Još jednom treba napomenuti da je razmer površina sračunat približno koristeći odnose površina konačnih figura, pa kod ove projekcije nije 1 jer nije sračunat po formuli za razmer površina. Kao što možemo primetiti, kod Sinusoidne projekcije imamo slučaj da presek između meridijana i paralela nije prav ugao.

tissot.sin <- st_transform(tissot.sf,"+proj=sinu +ellps=WGS84")
# površine u projekciji
tissot.sin$a_proj <- st_area(tissot.sin)
# odnos površina u projekciji i na elipsoidu
tissot.sin$p <- tissot.sin$a_proj/tissot.sin$a

pop = st_transform(pop, "+proj=sinu +ellps=WGS84")


ll = st_centroid(tissot.sin)
ll = st_coordinates(ll)

# Karta sa elipsama deformacija i razmerom površina

plot(st_geometry(pop),  graticule = TRUE, border = 'grey',  axes = TRUE)
plot(st_geometry(tissot.sin),  border = 'red', add = TRUE)
text(ll[, 1], ll[, 2], paste("p=",round(tissot.sin$p,1)) ,  cex = 0.8)

Slika 3.19: Elipse deformacija Sinusoidna projekcija.

3.6 Podela projekcija prema karakteru deformacija

Kao što smo videli, u opštem slučaju prilikom kartografskog preslikavanja deformišu se:

• DUŽINE,

• POVRŠINE,

• UGLOVI.

U pogledu deformacija projekcije se dele na:

• konformne (istougle, lokalno nema deformacije uglova, tj oblika),

• ekvivalente (istopovršinske) i

• uslovne (proizvoljne).

U grupu uslovnih spadaju i ekvidistantne, gde u jednom pravcu nema deformacije dužina. Najčešće u jednom glavnom pravcu; ili u pravcu meridijana; ili u pravcu paralele.

Projekcija uvek može pripadati samo jednoj grupi.

Konformne projekcije

Osim naziva konformne (istougle) projekcije, često se sreće i naziv ortomorfne projekcije, što znači da se ovim preslikavanjem obezbeđuje sličnost beskonačno malih figura ‒ tzv. lokalna sličnost.

U bilo kojoj tački istouglog preslikavanja indikatrisa je krug, tj. linearna razmera ostaje ista u svim pravcima oko zadate tačke, što se može izraziti odnosom:

\[ a=b=m=n=c \]

Ako se vratimo na formulu za deformaciju površina:

\[ p=\frac{\sqrt{E G} \sin (\theta) \mathrm{d} \lambda \mathrm{d} \varphi}{\mathrm{MN} \cos \varphi \mathrm{d} \lambda \mathrm{d} \varphi}=m n \sin (\theta) \]

Kod konformnih projekcija ugao između meridijana i parale je ={90}^o, linearna ramera je ista u svim pravcima pa imamo:

\[ p=a^2=b^2=m^2=n^2=c^2 \]

Istopovršinske (ekvivalentne projekcije)

Ekvivalentne (homalografske) projekcije očuvavaju površine figura i zadovoljavaju uslov:

\[ p=a b=m n \sin \theta=m n \cos \varepsilon=1 \]

Ako je razmer površine jednak jedinici u svakoj tački u ravni projekcije, onda i konačne površine ostaju očuvane i nema deformacije površina konačnih geometrijskih entiteta.

Uslovne projekcije

Kod ovih projekcija postoje istovremeno i deformacije uglova i deformacije površina. Ipak, iz veoma velike grupe uslovnih (proizvoljnih) projekcija izdvaja se grupa tzv. ekvidistantnih (istodužinskih) projekcija, kod kojih je linearan razmer uzduž jednog od glavnih pravaca konstantan i jednak jedinici, tj. ili a=1, ili b=1, kod ovih projekcija vrlo često se taj glavni pravac poklapa sa pravcem meridijana ili paralele.

3.7 Klasifikacija projekcija u odnosu na pomoćne projekcione površi

Zemljin elipsoid (ili lopta) ne samo da se preslikava direktno na ravan već se, ponekad, zamišljeno preslikava na bočne površi konusa ili cilindra, koje se zatim razvijaju u ravan, Slika 3.20.

1. Cilindrične,

2. Konusne i

3. Azimutne.

Slika 3.20: Klasifikacija projekcija u odnosu na pomoćne projekcione površi. Izvor: https://docs.qgis.org/3.28/en/_images/projection_families.png/

3.8 Klasifikacija projekcija u odnosu na položaj pomoćne projekcione površi

Položaj projekcionih površi iskazuje se, najočiglednije, položajem koji obrtne ose (cilindra ili konusa) zauzimaju u odnosu na obrtnu osu Zemlje, Slika 3.21. Kod azimutnih projekcija normala, na ravan projekcije i elipsoid (loptu) u dodirnoj tački. Prema tome, projekcije se dele na:

1. Normalne (prave),

2. Poprečne (transferzalne) i

3. Kose

Slika 3.21: Klasifikacija projekcija u odnosu na položaj (orijentaciju) pomoćne projekcione površi.

3.9 Klasifikacija projekcija u odnosu na zone

Jedinstvene projekcije

U ovom slučaju izrada karata za neko područje provodi se u jednom istom koordinatnom sistemu u ravni (npr. Merkatorova projekcija koju smo pominjali ranije (konformna), Sinusoidna (ekvivalentna), ili Robinsonova projekcija (ekvivalenta), Slika 3.22).

library(sf)
library(spData)

data(world)
granicaSRKS = st_union(world[173, 'geom' ], world[175,'geom'] )
world[173, ]$geom = granicaSRKS$geom
world = world[-175,]
world_r = st_transform(world, 
        crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84")
plot(world_r["lifeExp"], graticule = TRUE, 
     main = "Očekivani životni vek za rođene 2014.")

Slika 3.22: Robinsonova projekcija, jedinstvena globalna projekcija.

Višezone (višepojasne ) projekcije

Kod krupnorazmernih karata da bismo održali zadani nivo tačnosti, područje preslikavanja delimo na zone (duž meridijana), ili pojaseve (duž paralela). (Npr. UTM projekcija, prikaz 60 UTM zona, Slika 3.23.)

Slika 3.23: UTM zone. Izvor: https://snl.no/UTM

Višegrane projekcije

U ovom slučaju celokupno područje preslikavanja deli se linijama meridijana i paralela na odgovarajući broj sferoidnih ili sfernih trapeza. Svaki takav trapez preslikava se odvojeno u projekciji usvojenoj za celokupno područje preslikavanja. Svaki list karte je zaseban koordinatni sistem. (Modifikovana američka polikonusna projekcija, primer takve projekcije je i Nikolozi globular projekcija, koja je višegrana i deli ceo svet na 2 dela, Slika 3.24)

Slika 3.24: Nikolozi globular, višegrana projekcija sa dva prekida. Izvor: https://w.wiki/6$qZ

Prekinute projekcije

Ove se projekcije primenjuju za kartografisanje odvojenih delova Zemljine površi ‒ na primer, karte kontinenata ili vodenih površina, koje se spajaju po pravilu duž Ekvatora. Kartografska mreža prekinuta je duž odgovarajućih, unapred izabranih meridijana. Ponekad, štaviše, postoji ne samo jedan već nekoliko prekida mreže. (Gudeova (Goode homolosine) projekcija, Slika 3.25)

Slika 3.25: Gudeova projekcija sa prekidima. Izvor: https://en.wikipedia.org/wiki/Goode_homolosine_projection

Gimond, Manuel. 2022. Intro to GIS and Spatial Analysis. Colby College. Accessed June 2023. https://mgimond.github.io/Spatial/index.html.

Jovanović, Velibor. 1984. Matematička Kartografija. VGI.