U prethodnom delu knjige primeri su uglavnom bili primenjeni nad vektorskim podacima. U drugom delu knjige su bili primeri koji se odnose na prvi (direktni) i drugi (inverzni) zadatak matematičke kartografije zajedno sa pojedinim primerima koji uključuju i transformacije iz jednog koordinatnog sistema u drugi, ali nad vektorskim podacima. U ovom poglavlju će biti prikazana primena projekcija i koordinatnih transformacija nad rasterskim podacima, kao i kooridnatne transformacije iz jednog koordinatnog sistema u drugi u opštem slučaju.
10.1 Projekcije i transformacije rasterskih podataka
Kao što smo već videli kod vektorskih podataka, ako znamo da transformišemo jednu tačku, onda znamo i liniju koja predstavlja niz tačaka, isto tako i poligon koji je zatvorena linija i tako dalje.
Najpre kreirajmo jedan simulirani raster u gridu longitude i latitude, na WGS84 elipsoidu, koji je predstavljen dvodimenzionalno.
library(terra)# terra 1.7.3library(sf)# Warning: package 'sf' was built under R version 4.2.3# Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUEbbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan raster 30x30 stepenier <-rast(bbox, resolution=30 )er# class : SpatRaster # dimensions : 6, 12, 1 (nrow, ncol, nlyr)# resolution : 30, 30 (x, y)# extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)# coord. ref. : lon/lat WGS 84# vrednosti pikselaer$v=1:72plot(er)
Slika 10.1: Simulirani grid veličine 30x30 stepeni u gridu longitude i latitude na WGS84 elipsoidu.
Treba imati na umu da je prikazani grid u stvari bolje vizuelizovati u projekcijama kao što je Azimutna ortografska da bi se stekao utisak kako on izgleda na samom elipsoidu, Slika 10.2.
Slika 10.2: Izgled rastera 30x30 stepeni na Google Earth elipsoidu.
Sledeći primer slikovito pokazuje šta bi bilo kad bismo primenili takav pristup i kod rasterskih podataka. Kao što smo videli u prvom poglavlju, rasterske podatke možemo pretvoriti u vektorske pa bismo onda mogli transformisati vektorske podatke u željenu projekciju ili koordinatni referentni sistem, Slika 10.3. Dat je simulirani raster u gridu longituda latituda (\(30^o\)x\(30^o\)) koji je transformisan u poligonski vektorski podatak, a potom su poligoni piksela transformisani u Sinusoidnu projekciju..
Slika 10.3: Vektorizacija pa transformacija rastera u Sinusoidnu projekciju.
Očigledno prethodni primer ne funkcioniše. Originalno kod rasterskih podataka, pikseli su uređeni kao pravougaoni elementi (ili kvadratni), pri čemu svaki piksel ima istu geometriju. Transformisani raster je neophodno da ima pikselsku strukturu, gde svaki piksel takođe ima istu geometriju bilo da je ona kvadratna ili pravougaona. Dakle, neophodno je dobiti transformisani raster sa regularnom geometrijom piksela.
library(terra)library(sf)library(spData)# Warning: package 'spData' was built under R version 4.2.3bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan rasterer <-rast(bbox, resolution=30 )# vrednosti pikselaer$v=1:72data(world)granicaSRKS =st_union(world[173, 'geom' ], world[175,'geom'] )world[173, ]$geom = granicaSRKS$geomworld = world[-175,]# terra paket koristi svoj zapis za vektorske podatke lako je konvertovati # u sf objekat SpatVectorw_vect =vect(world)w_sin =project(w_vect,"+proj=sinu")# za terra paket transformacija se zadaje metodom project# rezolucija koju zadajemo u pravoglim koordinatame je oko 3000 km (1st oko 100km)er_sin =project(er, "+proj=sinu", res=30*100*1000, method="near")plot(er_sin, col=rev(terrain.colors(72)) )plot(w_sin, add=T)
Slika 10.4: Transformacija simuliranog rastera iz WGS84 grida u Sinusoidnu projekciju.
Kod transformacije rastera neophodno je da ponovo dobijemo rastersku strukturu, ali se pritom može pojaviti još jedan problem. Naime, pikseli se mogu preslikati i van područja koji obuhvata prikaz Zemlje (Wraparound Problem, Slika 10.4).
library(terra)library(sf)library(spData)bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan rasterer <-rast(bbox, resolution=30 )# vrednosti pikselaer$v=1:72# za terra paket transformacija se zadaje metodom project# rezolucija koju zadajemo u pravoglim koordinatame je oko 3000 km (1st oko 100km)# 30*100*1000er_sin=project(er, "+proj=sinu", res=30*100*1000, method="near")# zadržaćemo samo piksele koji se poklapaju sa vektorskim podacimapol =st_as_sf(as.polygons(er, dissolve=FALSE) )pol_sin =st_transform(pol, "+proj=sinu")er_sin = terra::mask(x = er_sin, mask =vect(pol_sin))g <-graticule(60, 30, crs="+proj=sinu")plot(er_sin, col=rev(terrain.colors(72)) )plot(as.polygons(er_sin, dissolve=FALSE) , border="gray", add=T)plot(w_sin, add=T)plot(g, col="blue",lty=3, retro=TRUE, add=T)
Slika 10.5: Transformacija simuliranog rastera iz WGS84 grida u Sinusoidnu projekciju.
Na prethodnim slikama jasno se vidi da pozicija i broj piksela nijesu isti u jednoj i drugoj projekciji, a da su vrednosti piksela preraspoređene tako da što vernije predstavljaju pozicije i vrednosti ulaznih piksela. Metode koje obezbeđuju način preraspoređivanja vrednosti piksela iz ulaznog rastera u izlazni raster nazivaju se resampling metode. Ovde će ukratko biti pomenute neke resampling metode. Kod svih gore navedenih primera korišćena je metoda najbližeg suseda project(er, "+proj=sinu", res=30*100*1000, method="near"), zadata argumentom method="near".
Resampling metode
Resampling metoda najbliži sused (Nearest Neighbor) je jednostavna metoda resamplinga. Ona uzima u obzir koji je najbliži centar piksela iz ulaznog rastera u odnosu na centre piksela izlaznog rastera, tako da nađe najbližeg suseda, kod takvog piksela dodeljuje se originalna vrednost ulaznog piksela. Ona se primenjuje kad je neophodno napraviti brzi resampling, kao i kod kategorijskih podataka, ili normalizovanih podataka.
Bilinearna interpolacija uzima u obzir 4 najbliže ćelije u pravcu glavnih osa rastera i računa vrednost u izlaznom rasteru kao srednju vrednost najbližih ćelija.
Kubna interpolacija koristi vrednosti 16 najbližih ćelija originalnog rastera da odredi vrednost izlazne ćelije, primenjujući polinomne funkcije trećeg reda. Koristi se za kontinuirane rastere i daje glatkiju površinu u poređenju sa bilinearnom interpolacijom, ali je računski zahtevnija.
Kubna interpolacija sa splajnom takođe koristi vrednosti 16 najbližih ćelija originalnog rastera da odredi vrednost izlazne ćelije, ali primenjuje kubne splajnove. Koristi se za kontinuirane rastere.
Uzmimo Digitalni model visina (DEM) u okolini Zlatibora, podaci se mogu preuzeti iz paketageodata .
library(terra)library(sf)library(geodata)# Warning: package 'geodata' was built under R version 4.2.3# alt = elevation_30s('SRB', path=tempdir() )# terra::res(alt)# 0.008333333 rezolucija oko 1km rezolucija # 0.008333333*pi/180*6377 km je 0.93 kmalt =rast("/vsicurl/https://zenodo.org/record/6211553/files/copernicus_DSM_30arcsec_mood_bbox_epsg4326.tif") # https://zenodo.org/record/6211553names(alt) ="h"zlatibor =crop(alt, ext(19.65,19.75, 43.60, 43.75)) # meters * 1000 (scaled to Integer; example: value 23220 = 23.220 m a.s.l.)values(zlatibor) =round(values(zlatibor)/1000)pol =st_as_sf(as.polygons(zlatibor, dissolve=FALSE) )plot(zlatibor, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( pol$h) , cex=0.4 )
Slika 10.6: Prikaz DEM-a u okolini Zlatibora.
Transfromisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Najbližeg suseda gde je rezolucija grida u pravouglim koordinatama 1.000 m, Slika 10.7.
Slika 10.7: Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Najbliži sused.
Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Bilinearne interpolacije gde je rezolucija grida u pravouglim koordinatama 1.000 m, Slika 10.8.
Slika 10.8: Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Bilinearna interpolacija.
Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Kubne interpolacije, gde je rezolucija grida u pravouglim koordinatama 1.000 m, Slika 10.9.
Slika 10.9: Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Kubna interpolacija.
Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Kubne splajn interpolacije, gde je rezolucija grida u pravouglim koordinatama 1.000 m, Slika 10.10.
Slika 10.10: Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Kubna interpolacija sa splajnom.
Na prethodnim kartama možemo uočiti beli piksel koji ima najveću vrednost i rezultat resampling metoda koje preraspoređuju vrednosti piksela u novu geometriju. Što je veće glačanje, vrednost piksela sa najvećom visinom je manja, jer u proračun ulaze okolni pikseli, osim kod metode Najbižeg suseda, videti: Slika 10.7, Slika 10.8, Slika 10.9, Slika 10.10.
Kao što je već napomenuto u prvom poglavlju, GDAL funkcionalnosti se često koriste za razne vrste rasterskih funkcionalnosti i mogu se koristiti i direktno iz komandne linije ili u okviru paketa implementiranih u: R-u, Pajtonu ili okviru klasičnih GIS softvera, npr. QGIS-a. Primer korišćenja gdalwarp funkcije ilustrovan je na podacima Natural Earth koji predstavljaju kontinentalne mase i vodene površine, stilizacija je integrisana u TIFF fajlu.
Parametar -t_srs koristi se za zadavanje koordinatnog sistema izlaznog fajla, -r se koristi za zadavanje metoda resamplinga, ovde je zadata kubna konvolucija, -co parametar za izbor kompresije, -dstalpha definiše da se pikseli koji nemaju vrednost prikažu transparentno. Ulazni fajl je NE1_50M_SR_W.tif, prikazan je na slici (Slika 10.11). Kod rasterskih funkcija, za razliku od vektorskih, izlazni fajl je uvek poslednji parametar funkcije, NE1_50M_SR_W_sin.tif, Slika 10.12).
Slika 10.11: Prikaz kontinentalnih masa u WGS84 koordinatnom sistemu NE150MSRW.tif.
Slika 10.12: Prikaz kontinentalnih masa u Sinusoidnoj projekciji NE150MSRWsin.tif.
10.2 Helmertova i gridna koordinatna transformacija
Ranije smo prolazili kroz primere transformacije iz jednog koordinatnog sistema u drugi, uglavnom se u pozadini primenjivala Helmertova 7-parametarska transformacija. Može se desiti da za neke podatke imamo i datum definisan samo sa 3 parametra translacije, a ne svih 7 koji obuhvataju i translacije, i rotacije i razmeru. Sledeći primer ilustruje razlike u računanju iste tačke ako koristimo 2 različita datuma.
Elipsoidne koordinate spomenika Nikoli Tesli ispred Građevinskog fakulteta u Beogradu (\(\varphi=44.8053183\) i \(\lambda=20.4763281\)) iz koordinatnog sistema WGS84 transformisati u Gaus-Krigerovu projekciju za datum koji je definisan samo sa tri translacije +towgs84=682,-203,480,0,0,0,0 i u gedetski datum zadat sa svih 7 parametara +towgs84=574.027,170.175,401.545,4.88786,-0.66524,-13.24673,6.8893.
library(sf)# globalna opcija za broj cifara koji se ispisuju u konzoli options(digits =12)# Kreirati prostorni lejer sa pozicijom zgrade Građevinskog fakulteta # prvi nivo je je sfg objekat koji je samo geometrijski primitivgrf =st_point( c(20.4763281, 44.8053183) ) # lon , lat# potom pravimo kolekciju geometrije sfc u odgovarajućoj projekcijipts_sfc =st_sfc(grf, crs=4326) #epsg 4326 je WGS84 long lat# GK Srb koord. ref. sis. 3 translacijegk3 ="+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs +type=crs"# EPSG 3909# GK Srb koord. ref. sis. 7 parametragk7 ="+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=574.027,170.175,401.545,4.88786,-0.66524,-13.24673,6.8893 +units=m"# transformacija u UTMpts_sfc_gk3 =st_transform(pts_sfc, crs =gk3)pts_sfc_gk7 =st_transform(pts_sfc, crs =gk7)# matrica koordinata iz prostornih podataka 3 translacijest_coordinates(pts_sfc_gk3)# X Y# [1,] 7459008.35955 4962434.20491st_coordinates(pts_sfc_gk7)# X Y# [1,] 7459006.31365 4962438.72797# NAPOMENA u R sf paketu koordinate imaju oznaku X,Y kao u matematici# REDOSLED KOORDINATA JE E,N (pravac istoka, pravac severa)
Iz navedenog primera je jasno da parametri geodetskog datuma znatno utiču na tačnost koordinatnih transformacija.
Gridna transformacija je implementirana u novijim verzijama biblioteke PROJ koja je implementirana u R-u. Izgled datumskog grida za horizontalnu transformaciju dat je na sledećoj slici, Slika 10.13.
OSTN15 je definitivna transformacija OSGB36/ETRS89, https://github.com/OSGeo/PROJ-data/tree/master/uk_os. OSTN15 u kombinaciji sa ETRS89 koordinatama OS Net stanica, umesto stare triangulacione mreže, definišu nacionalnu mrežu. To znači da će, na primer, koordinate nacionalne mreže postojeće OSGB36 tačke, rekalibrisane korišćenjem GNSS-a iz OS Net i OSTN15, biti ispravne.
library(sf)library(terra)grid <-rast("/vsicurl/https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif")# u sekundamaplot(grid[[1:2]])
Slika 10.13: Horizontalni datumski grid za Veliku Britaniju za transformaciju iz OSGB36 u ETRS89, vrednosti su date lučnim sekundama.
Sledeći primer demonstrira koordinatnu transformaciju koja koristi gridni datum.
Primer vertikalnog datumskog grida dat je na sledećoj karti, Slika 10.14. Vertikalna gridna transformacija za Geoid model OSGM15, ODN visina (EPSG:5701). Koristi se u transformaciji sa OSGM15 ortometrijskih visina na ETRS89 elipsoidne visine. Raster praktično predstavlja undulacije geoida, odnosno geoidne visine.
library(sf)library(terra)vgrid <-rast("/vsicurl/https://cdn.proj.org/uk_os_OSGM15_GB.tif")# u metrimaplot(vgrid)
Slika 10.14: Vertikalni datumski grid za Veliku Britaniju u metrima, undulacije geoida.
Primer zadavanja vertikalne transformacije u proj notaciji:
Afina transformacija je transformacija između afinih prostora koja preslikava tačke u tačke, prave u prave i ravni u ravni. Takođe, kod afinih preslikavanja par paralelnih pravih ostaje paralelan posle primene transformacije, ali afina transformacija ne mora nužno da sačuva uglove između pravih ili rastojanja između tačaka, mada čuva razmeru kolinearnih tačaka. Afine transformacije uključuju, između ostalog, translacije (pomeranje), skaliranje (urazmeravanje) i rotaciju u prostoru. Pored toga, moguće je koristiti bilo koju kombinaciju, nrp.: samo translacija i skaliranje, translacija i rotiranje, samo translacija, samo rotiranje ili sve zajedno. Često se afina transformacija primenjuje kao pomoćno sredstvo u kartografiji, geoinformatici i GIS-u, npr: kod pozicioniranja toponima često se koristi samo translacija. Afina transformacija se koristi i kod transformacije podataka iz starih skeniranih karata (npr. iz Stereografske projekcije 1:2880 u DKS na području Vojvodine). Primenjuje se u fotogrametriji. Afina transformacija se primenjuje i kad imamo pogrešnu geometriju prostornih podataka koju treba ispraviti itd. Paket sf implementira afinu transformaciju za objekte klasa sfg i sfc.
Gde su \(x,y\) koordinate iz jednog koordinatnog sistema, \(x',y'\) koordinate iz drugog koordinatnog sistema (transformisanog). Malim slovima \(a\) i \(b\) označeni su parametri transformacije.
Prethodna formula se može izraziti i u obliku:
\[
\left[\begin{array}{c}
x^{\prime} \\
y^{\prime} \\
1
\end{array}\right]=\left[\begin{array}{ccc}
a_{11} & a_{12} & b_1 \\
a_{21} & a_{22} & b_2 \\
0 & 0 & 1
\end{array}\right]\left[\begin{array}{l}
x \\
y \\
1
\end{array}\right]
\] Dakle, kod afine transformacije u opštem slučaju imamo 6 parametara transformacije (označimo ih na sledeći način \(a, b, c\) za koordinatu u pravcu istoka i \(d,f,g\) koji se odnose na koordinatu u pravcu severa u geodeziji \(y,x\)), ako želimo da postavimo model za metod najmanjih kvadrata, kakav je uobičajen u geodeziji, kojim bismo rešili vrednosti parametara kad imamo grupe zajedničkih tačaka, imali bismo:
\[
\begin{aligned}
& a \cdot y+b \cdot x+c=y'+v_y \\
& d \cdot y+f \cdot x+g=x'+v_x
\end{aligned}
\] Recimo da imamo 4 identične tačke, važile bi jednačine:
U geodeziji je uobičajeno da afinu transformaciju definišemo na sledeći način kad imamo definisan jedan parametar razmere (skaliranja) koji je isti za obe koordinatne ose, a koordinate su zarotirane za ugao \(\theta\):
\[
\left[\begin{array}{c}
X \\
Y \\
Z \\
T
\end{array}\right]^{\text {dest }}=\left[\begin{array}{c}
\text { xoff } \\
\text { yoff } \\
\text { zoff } \\
\text { toff }
\end{array}\right]+\left[\begin{array}{cccc}
s 11 & s 12 & s 13 & 0 \\
s 21 & s 22 & s 23 & 0 \\
s 31 & s 32 & s 33 & 0 \\
0 & 0 & 0 & \text { tscale }
\end{array}\right]\left[\begin{array}{c}
X \\
Y \\
Z \\
T
\end{array}\right]^{\text {source }}
\]
Iz dokumentacije sf::vignette("sf3")(Pebesma 2018) primer ilustruje afinu transformaciju. Sve države Severne Karoline su zarotirane za \(90^o\) oko centroida svakog poligona i skalirane faktorom 0.75, Slika 10.15.
# Rad sa rasterima i koordinatne transformacije {#sec-trans}U prethodnom delu knjige primeri su uglavnom bili primenjeni nad vektorskim podacima. U drugom delu knjige su bili primeri koji se odnose na prvi (direktni) i drugi (inverzni) zadatak matematičke kartografije zajedno sa pojedinim primerima koji uključuju i transformacije iz jednog koordinatnog sistema u drugi, ali nad vektorskim podacima. U ovom poglavlju će biti prikazana primena projekcija i koordinatnih transformacija nad rasterskim podacima, kao i kooridnatne transformacije iz jednog koordinatnog sistema u drugi u opštem slučaju.## Projekcije i transformacije rasterskih podatakaKao što smo već videli kod vektorskih podataka, ako znamo da transformišemo jednu tačku, onda znamo i liniju koja predstavlja niz tačaka, isto tako i poligon koji je zatvorena linija i tako dalje.Najpre kreirajmo jedan simulirani raster u gridu longitude i latitude, na WGS84 elipsoidu, koji je predstavljen dvodimenzionalno.```{r fig-r0, echo=!knitr::is_latex_output() }#| out.width: 60%#| fig.cap: "Simulirani grid veličine 30x30 stepeni u gridu longitude i latitude na WGS84 elipsoidu."#| code-fold: falselibrary(terra)library(sf)bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan raster 30x30 stepenier <-rast(bbox, resolution=30 )er# vrednosti pikselaer$v=1:72plot(er)```Treba imati na umu da je prikazani grid u stvari bolje vizuelizovati u projekcijama kao što je Azimutna ortografska da bi se stekao utisak kako on izgleda na samom elipsoidu, @fig-rgee.```{r fig-rgee, out.width='100%' , fig.align = "center", out.extra = '', fig.cap='Izgled rastera 30x30 stepeni na Google Earth elipsoidu.', echo=FALSE}knitr::include_graphics('images/geerast.jpg')```Sledeći primer slikovito pokazuje šta bi bilo kad bismo primenili takav pristup i kod rasterskih podataka. Kao što smo videli u prvom poglavlju, rasterske podatke možemo pretvoriti u vektorske pa bismo onda mogli transformisati vektorske podatke u željenu projekciju ili koordinatni referentni sistem, @fig-r1. Dat je simulirani raster u gridu longituda latituda ($30^o$x$30^o$) koji je transformisan u poligonski vektorski podatak, a potom su poligoni piksela transformisani u Sinusoidnu projekciju..```{r fig-r1, echo=!knitr::is_latex_output() }#| fig.cap: "Vektorizacija pa transformacija rastera u Sinusoidnu projekciju."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymxn, ymax)# prazan raster 30x30 stepenier <-rast(bbox, resolution=30 ) # vrednosti pikselaer$v=1:72pol =st_as_sf(as.polygons(er, dissolve=FALSE) )layout(matrix(1:2, 1, 2), c(1,1))par(mar =rep(0, 4))plot(st_geometry(pol), col=rev(terrain.colors(72)) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste(pol$v) , cex=0.6 )pol_sin =st_transform(pol, "+proj=sinu")plot(st_geometry(pol_sin), col=rev(terrain.colors(72)) )cc =st_coordinates(st_centroid(st_geometry(pol_sin)))text(cc[, 1], cc[, 2], labels =paste(pol$v) , cex=0.6 )```Očigledno prethodni primer ne funkcioniše. Originalno kod rasterskih podataka, pikseli su uređeni kao pravougaoni elementi (ili kvadratni), pri čemu svaki piksel ima istu geometriju. Transformisani raster je neophodno da ima pikselsku strukturu, gde svaki piksel takođe ima istu geometriju bilo da je ona kvadratna ili pravougaona. Dakle, neophodno je dobiti transformisani raster sa regularnom geometrijom piksela.```{r fig-r2, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija simuliranog rastera iz WGS84 grida u Sinusoidnu projekciju."#| out.width: 90%#| code-fold: falselibrary(terra)library(sf)library(spData)bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan rasterer <-rast(bbox, resolution=30 )# vrednosti pikselaer$v=1:72data(world)granicaSRKS =st_union(world[173, 'geom' ], world[175,'geom'] )world[173, ]$geom = granicaSRKS$geomworld = world[-175,]# terra paket koristi svoj zapis za vektorske podatke lako je konvertovati # u sf objekat SpatVectorw_vect =vect(world)w_sin =project(w_vect,"+proj=sinu")# za terra paket transformacija se zadaje metodom project# rezolucija koju zadajemo u pravoglim koordinatame je oko 3000 km (1st oko 100km)er_sin =project(er, "+proj=sinu", res=30*100*1000, method="near")plot(er_sin, col=rev(terrain.colors(72)) )plot(w_sin, add=T)```Kod transformacije rastera neophodno je da ponovo dobijemo rastersku strukturu, ali se pritom može pojaviti još jedan problem. Naime, pikseli se mogu preslikati i van područja koji obuhvata prikaz Zemlje (*Wraparound Problem*, @fig-r2).```{r fig-r3, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija simuliranog rastera iz WGS84 grida u Sinusoidnu projekciju."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)library(spData)bbox=ext(c(-180, 180, -90, 90) )#(xmin, xmax, ymin, ymax)# prazan rasterer <-rast(bbox, resolution=30 )# vrednosti pikselaer$v=1:72# za terra paket transformacija se zadaje metodom project# rezolucija koju zadajemo u pravoglim koordinatame je oko 3000 km (1st oko 100km)# 30*100*1000er_sin=project(er, "+proj=sinu", res=30*100*1000, method="near")# zadržaćemo samo piksele koji se poklapaju sa vektorskim podacimapol =st_as_sf(as.polygons(er, dissolve=FALSE) )pol_sin =st_transform(pol, "+proj=sinu")er_sin = terra::mask(x = er_sin, mask =vect(pol_sin))g <-graticule(60, 30, crs="+proj=sinu")plot(er_sin, col=rev(terrain.colors(72)) )plot(as.polygons(er_sin, dissolve=FALSE) , border="gray", add=T)plot(w_sin, add=T)plot(g, col="blue",lty=3, retro=TRUE, add=T)```Na prethodnim slikama jasno se vidi da pozicija i broj piksela nijesu isti u jednoj i drugoj projekciji, a da su vrednosti piksela preraspoređene tako da što vernije predstavljaju pozicije i vrednosti ulaznih piksela. Metode koje obezbeđuju način preraspoređivanja vrednosti piksela iz ulaznog rastera u izlazni raster nazivaju se **resampling** metode. Ovde će ukratko biti pomenute neke resampling metode. Kod svih gore navedenih primera korišćena je metoda najbližeg suseda `project(er, "+proj=sinu", res=30*100*1000, method="near")`, zadata argumentom ` method="near"`.::: callout-note#### Resampling metode**Resampling metoda najbliži sused** (*Nearest Neighbor*) je jednostavna metoda resamplinga. Ona uzima u obzir koji je najbliži centar piksela iz ulaznog rastera u odnosu na centre piksela izlaznog rastera, tako da nađe najbližeg suseda, kod takvog piksela dodeljuje se originalna vrednost ulaznog piksela. Ona se primenjuje kad je neophodno napraviti brzi resampling, kao i kod kategorijskih podataka, ili normalizovanih podataka.**Bilinearna interpolacija** uzima u obzir 4 najbliže ćelije u pravcu glavnih osa rastera i računa vrednost u izlaznom rasteru kao srednju vrednost najbližih ćelija.**Kubna interpolacija** koristi vrednosti 16 najbližih ćelija originalnog rastera da odredi vrednost izlazne ćelije, primenjujući polinomne funkcije trećeg reda. Koristi se za kontinuirane rastere i daje glatkiju površinu u poređenju sa bilinearnom interpolacijom, ali je računski zahtevnija.**Kubna interpolacija sa splajnom** takođe koristi vrednosti 16 najbližih ćelija originalnog rastera da odredi vrednost izlazne ćelije, ali primenjuje kubne splajnove. Koristi se za kontinuirane rastere.:::Uzmimo Digitalni model visina (DEM) u okolini Zlatibora, podaci se mogu preuzeti iz paketa`geodata` .```{r fig-dem, echo=!knitr::is_latex_output()}#| fig.cap: "Prikaz DEM-a u okolini Zlatibora."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)library(geodata)# alt = elevation_30s('SRB', path=tempdir() )# terra::res(alt)# 0.008333333 rezolucija oko 1km rezolucija # 0.008333333*pi/180*6377 km je 0.93 kmalt =rast("/vsicurl/https://zenodo.org/record/6211553/files/copernicus_DSM_30arcsec_mood_bbox_epsg4326.tif") # https://zenodo.org/record/6211553names(alt) ="h"zlatibor =crop(alt, ext(19.65,19.75, 43.60, 43.75)) # meters * 1000 (scaled to Integer; example: value 23220 = 23.220 m a.s.l.)values(zlatibor) =round(values(zlatibor)/1000)pol =st_as_sf(as.polygons(zlatibor, dissolve=FALSE) )plot(zlatibor, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( pol$h) , cex=0.4 )```Transfromisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Najbližeg suseda gde je rezolucija grida u pravouglim koordinatama 1.000 m, @fig-resn.```{r fig-resn, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Najbliži sused."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)utm_srb ="+proj=utm +zone=34 +ellps=GRS80 +towgs84=0.26901,0.18246,0.06872,-0.01017,0.00893,-0.01172,0.04 +units=m"# transformacija u UTM 1000 m rez metod neardem_n =project(zlatibor, utm_srb, res=1000, method ="near")pol =st_as_sf(as.polygons(dem_n, dissolve=FALSE) )plot(dem_n, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( pol$h) , cex=0.4 )```Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Bilinearne interpolacije gde je rezolucija grida u pravouglim koordinatama 1.000 m, @fig-resb.```{r fig-resb, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Bilinearna interpolacija."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)utm_srb ="+proj=utm +zone=34 +ellps=GRS80 +towgs84=0.26901,0.18246,0.06872,-0.01017,0.00893,-0.01172,0.04 +units=m"# transformacija u UTM 1000 m rez metod bilinear dem_b =project(zlatibor, utm_srb, res=1000, method ="bilinear")pol =st_as_sf(as.polygons(dem_b, dissolve=FALSE) )plot(dem_b, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( round(pol$h,1) ) , cex=0.4 )```Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Kubne interpolacije, gde je rezolucija grida u pravouglim koordinatama 1.000 m, @fig-resc.```{r fig-resc, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Kubna interpolacija."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)utm_srb ="+proj=utm +zone=34 +ellps=GRS80 +towgs84=0.26901,0.18246,0.06872,-0.01017,0.00893,-0.01172,0.04 +units=m"# transformacija u UTM 1000 m rez metod cubicdem_c =project(zlatibor, utm_srb, res=1000, method ="cubic")pol =st_as_sf(as.polygons(dem_c, dissolve=FALSE) )plot(dem_c, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( round(pol$h,1) ) , cex=0.4 )```Transformisati DEM iz WGS84 koordinatnog sistema u UTM projekciju za Srbiju koristeći resampling metod Kubne splajn interpolacije, gde je rezolucija grida u pravouglim koordinatama 1.000 m, @fig-rescs.```{r fig-rescs, echo=!knitr::is_latex_output() }#| fig.cap: "Transformacija DEM rastera iz WGS84 grida u UTM projekciju, metod Kubna interpolacija sa splajnom."#| out.width: 100%#| code-fold: falselibrary(terra)library(sf)utm_srb ="+proj=utm +zone=34 +ellps=GRS80 +towgs84=0.26901,0.18246,0.06872,-0.01017,0.00893,-0.01172,0.04 +units=m"# transformacija u UTM 1000 m rez metod cubicsplinedem_cs =project(zlatibor, utm_srb, res=1000, method ="cubicspline")pol =st_as_sf(as.polygons(dem_cs, dissolve=FALSE) )plot(dem_cs, col=terrain.colors(160,alpha=0.5) )cc =st_coordinates(st_centroid(st_geometry(pol)))text(cc[, 1], cc[, 2], labels =paste( round(pol$h,1) ) , cex=0.4 )```Na prethodnim kartama možemo uočiti beli piksel koji ima najveću vrednost i rezultat resampling metoda koje preraspoređuju vrednosti piksela u novu geometriju. Što je veće glačanje, vrednost piksela sa najvećom visinom je manja, jer u proračun ulaze okolni pikseli, osim kod metode Najbižeg suseda, videti: @fig-resn, @fig-resb, @fig-resc, @fig-rescs.Kao što je već napomenuto u prvom poglavlju, GDAL funkcionalnosti se često koriste za razne vrste rasterskih funkcionalnosti i mogu se koristiti i direktno iz komandne linije ili u okviru paketa implementiranih u: R-u, Pajtonu ili okviru klasičnih GIS softvera, npr. QGIS-a. Primer korišćenja [`gdalwarp`](http://www.gdal.org/gdalwarp.html) funkcije ilustrovan je na podacima [Natural Earth](http//www.naturalearthdata.com/download/50m/raster/NE1_50M_SR_W.zip) koji predstavljaju kontinentalne mase i vodene površine, stilizacija je integrisana u TIFF fajlu.<!-- \small -->```{bash, out.width='90%', eval=FALSE}gdalwarp-t_srs'+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'-r cubic -co COMPRESS=LZW -dstalpha NE1_50M_SR_W.tif NE1_50M_SR_W_sin.tif```\normalsizeParametar `-t_srs` koristi se za zadavanje koordinatnog sistema izlaznog fajla, `-r` se koristi za zadavanje metoda resamplinga, ovde je zadata kubna konvolucija, `-co` parametar za izbor kompresije, `-dstalpha` definiše da se pikseli koji nemaju vrednost prikažu transparentno. Ulazni fajl je `NE1_50M_SR_W.tif`, prikazan je na slici (@fig-ne1). Kod rasterskih funkcija, za razliku od vektorskih, izlazni fajl je uvek poslednji parametar funkcije, `NE1_50M_SR_W_sin.tif`, @fig-ne2).```{r fig-ne1, out.width='60%' , fig.align = "center" , out.extra = '', fig.cap='Prikaz kontinentalnih masa u WGS84 koordinatnom sistemu NE150MSRW.tif.', echo=FALSE}knitr::include_graphics('images/_media1/image36.png')``````{r fig-ne2, out.width='60%' , fig.align = "center" , out.extra = '', fig.cap='Prikaz kontinentalnih masa u Sinusoidnoj projekciji NE150MSRWsin.tif.', echo=FALSE}knitr::include_graphics('images/_media1/image37.png')```## Helmertova i gridna koordinatna transformacijaRanije smo prolazili kroz primere transformacije iz jednog koordinatnog sistema u drugi, uglavnom se u pozadini primenjivala Helmertova 7-parametarska transformacija. Može se desiti da za neke podatke imamo i datum definisan samo sa 3 parametra translacije, a ne svih 7 koji obuhvataju i translacije, i rotacije i razmeru. Sledeći primer ilustruje razlike u računanju iste tačke ako koristimo 2 različita datuma.Elipsoidne koordinate spomenika Nikoli Tesli ispred Građevinskog fakulteta u Beogradu ($\varphi=44.8053183$ i $\lambda=20.4763281$) iz koordinatnog sistema WGS84 transformisati u Gaus-Krigerovu projekciju za datum koji je definisan samo sa tri translacije `+towgs84=682,-203,480,0,0,0,0` i u gedetski datum zadat sa svih 7 parametara `+towgs84=574.027,170.175,401.545,4.88786,-0.66524,-13.24673,6.8893`. ```{r sfutm1, message=FALSE, warning=FALSE}library(sf)# globalna opcija za broj cifara koji se ispisuju u konzoli options(digits =12)# Kreirati prostorni lejer sa pozicijom zgrade Građevinskog fakulteta # prvi nivo je je sfg objekat koji je samo geometrijski primitivgrf =st_point( c(20.4763281, 44.8053183) ) # lon , lat# potom pravimo kolekciju geometrije sfc u odgovarajućoj projekcijipts_sfc =st_sfc(grf, crs=4326) #epsg 4326 je WGS84 long lat# GK Srb koord. ref. sis. 3 translacijegk3 ="+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=682,-203,480,0,0,0,0 +units=m +no_defs +type=crs"# EPSG 3909# GK Srb koord. ref. sis. 7 parametragk7 ="+proj=tmerc +lat_0=0 +lon_0=21 +k=0.9999 +x_0=7500000 +y_0=0 +ellps=bessel +towgs84=574.027,170.175,401.545,4.88786,-0.66524,-13.24673,6.8893 +units=m"# transformacija u UTMpts_sfc_gk3 =st_transform(pts_sfc, crs =gk3)pts_sfc_gk7 =st_transform(pts_sfc, crs =gk7)# matrica koordinata iz prostornih podataka 3 translacijest_coordinates(pts_sfc_gk3)st_coordinates(pts_sfc_gk7)# NAPOMENA u R sf paketu koordinate imaju oznaku X,Y kao u matematici# REDOSLED KOORDINATA JE E,N (pravac istoka, pravac severa)```Iz navedenog primera je jasno da parametri geodetskog datuma znatno utiču na tačnost koordinatnih transformacija.Gridna transformacija je implementirana u novijim verzijama biblioteke **PROJ** koja je implementirana u R-u. Izgled datumskog grida za horizontalnu transformaciju dat je na sledećoj slici, @fig-datumg.OSTN15 je definitivna transformacija OSGB36/ETRS89, <https://github.com/OSGeo/PROJ-data/tree/master/uk_os>. OSTN15 u kombinaciji sa ETRS89 koordinatama OS Net stanica, umesto stare triangulacione mreže, definišu nacionalnu mrežu. To znači da će, na primer, koordinate nacionalne mreže postojeće OSGB36 tačke, rekalibrisane korišćenjem GNSS-a iz OS Net i OSTN15, biti ispravne.```{r fig-datumg, message=FALSE, warning=FALSE}#| fig.cap: "Horizontalni datumski grid za Veliku Britaniju za transformaciju iz OSGB36 u ETRS89, vrednosti su date lučnim sekundama."#| out.width: 100%#| code-fold: falselibrary(sf)library(terra)grid <-rast("/vsicurl/https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif")# u sekundamaplot(grid[[1:2]])```Sledeći primer demonstrira koordinatnu transformaciju koja koristi gridni datum.```{r datumg2, message=FALSE, warning=FALSE}options(digits =12)library(sf)sf_proj_network(TRUE)b_pump <-st_read(system.file("gpkg/b_pump.gpkg", package="sf"))sf_proj_network()# EPSG 27700# +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 # +y_0=-100000 +ellps=airy +nadgrids=OSTN15_NTv2_OSGBtoETRS.gsb # +units=m +no_defs +type=crsbp <-sf_proj_pipelines("EPSG:27700", "EPSG:32630")bp# Gridna transformacijautm <-st_transform(b_pump, "EPSG:32630")st_coordinates(utm)# sa helmertovom transformacijomuk_7p="+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +nodefs"b_pump_7p = b_pumpst_crs(b_pump_7p ) = uk_7putm2 <-st_transform(b_pump_7p , "EPSG:32630")st_coordinates(utm2)```Primer vertikalnog datumskog grida dat je na sledećoj karti, @fig-vertgrid. Vertikalna gridna transformacija za Geoid model OSGM15, ODN visina (EPSG:5701). Koristi se u transformaciji sa OSGM15 ortometrijskih visina na ETRS89 elipsoidne visine. Raster praktično predstavlja undulacije geoida, odnosno geoidne visine.```{r fig-vertgrid, message=FALSE, warning=FALSE}#| fig.cap: "Vertikalni datumski grid za Veliku Britaniju u metrima, undulacije geoida."#| out.width: 100%#| code-fold: falselibrary(sf)library(terra)vgrid <-rast("/vsicurl/https://cdn.proj.org/uk_os_OSGM15_GB.tif")# u metrimaplot(vgrid)```Primer zadavanja vertikalne transformacije u `proj` notaciji:`+proj=vgridshift +grids=uk_os_OSGM15_GB.tif +multiplier=1`Ovako zadata transformacija primenjuje se na sledeći način za gridne (rasterske podatke):$$Z_{\text {target }}=Z_{\text {source }}+\text { multiplier } \times \text { gridvalue }$$## Afina transformacijaAfina transformacija je transformacija između afinih prostora koja preslikava tačke u tačke, prave u prave i ravni u ravni. Takođe, kod afinih preslikavanja **par paralelnih pravih ostaje paralelan posle primene transformacije**, ali afina transformacija ne mora nužno da sačuva uglove između pravih ili rastojanja između tačaka, mada čuva razmeru kolinearnih tačaka. Afine transformacije uključuju, između ostalog, translacije (pomeranje), skaliranje (urazmeravanje) i rotaciju u prostoru. Pored toga, moguće je koristiti bilo koju kombinaciju, nrp.: samo translacija i skaliranje, translacija i rotiranje, samo translacija, samo rotiranje ili sve zajedno. Često se afina transformacija primenjuje kao pomoćno sredstvo u kartografiji, geoinformatici i GIS-u, npr: kod pozicioniranja toponima često se koristi samo translacija. Afina transformacija se koristi i kod transformacije podataka iz starih skeniranih karata (npr. iz Stereografske projekcije 1:2880 u DKS na području Vojvodine). Primenjuje se u fotogrametriji. Afina transformacija se primenjuje i kad imamo pogrešnu geometriju prostornih podataka koju treba ispraviti itd. Paket sf implementira afinu transformaciju za objekte klasa `sfg` i `sfc`.Afina transformacija u ravni ima oblik:$$\left[\begin{array}{l}x^{\prime} \\y^{\prime}\end{array}\right]=\left[\begin{array}{ll}a_{11} & a_{12} \\a_{21} & a_{22}\end{array}\right]\left[\begin{array}{l}x \\y\end{array}\right]+\left[\begin{array}{l}b_1 \\b_2\end{array}\right]$$Gde su $x,y$ koordinate iz jednog koordinatnog sistema, $x',y'$ koordinate iz drugog koordinatnog sistema (transformisanog). Malim slovima $a$ i $b$ označeni su parametri transformacije.Prethodna formula se može izraziti i u obliku:$$\left[\begin{array}{c}x^{\prime} \\y^{\prime} \\1\end{array}\right]=\left[\begin{array}{ccc}a_{11} & a_{12} & b_1 \\a_{21} & a_{22} & b_2 \\0 & 0 & 1\end{array}\right]\left[\begin{array}{l}x \\y \\1\end{array}\right]$$Dakle, kod afine transformacije u opštem slučaju imamo 6 parametara transformacije (označimo ih na sledeći način $a, b, c$ za koordinatu u pravcu istoka i $d,f,g$ koji se odnose na koordinatu u pravcu severa u geodeziji $y,x$), ako želimo da postavimo model za metod najmanjih kvadrata, kakav je uobičajen u geodeziji, kojim bismo rešili vrednosti parametara kad imamo grupe zajedničkih tačaka, imali bismo:$$\begin{aligned}& a \cdot y+b \cdot x+c=y'+v_y \\& d \cdot y+f \cdot x+g=x'+v_x\end{aligned}$$Recimo da imamo 4 identične tačke, važile bi jednačine:<!-- $$ --><!-- \left.\begin{array}{l} --><!-- a \cdot y_1+b \cdot x_1+c=y_1^{\prime}+v_{y_1} \\ --><!-- d \cdot y_1+\mathrm{f} \cdot x_1^{\prime}+g=x_1^{\prime}+v_{x_1} --><!-- \end{array} \right\} \quad \text(T_1) --><!-- $$ --><!-- $$ --><!-- \left.\begin{array}{l}a \cdot y_2+b \cdot x_2+c=y'_2+v_{y_2} \\ --><!-- d \cdot y_2+\mathrm{f} \cdot x_2+g=x'_2+v_{x_2}\end{array}\right\} \quad \text(T_2) --><!-- $$ --><!-- $$ --><!-- \left.\begin{array}{l}a \cdot y_3+b \cdot x_3+c=y'_3+v_{y_3} \\ --><!-- d \cdot y_3+\mathrm{f} \cdot x_3+g=x'_3+v_{x_3}\end{array}\right\} \quad \text(T_3) --><!-- $$ --><!-- $$ --><!-- \left.\begin{array}{l}a \cdot y_4+b \cdot x_4+c=y'_4+v_{y_4} \\ --><!-- d \cdot y_4+\mathrm{f} \cdot x_4+g=x'_4+v_{x_4}\end{array}\right\} \quad \text(T_4) --><!-- $$ -->$$\begin{aligned}&a \cdot y_1+b \cdot x_1+c=y_1^{\prime}+v_{y_1} \\&d \cdot y_1+\mathrm{f} \cdot x_1^{\prime}+g=x_1^{\prime}+v_{x_1} \\&a \cdot y_2+b \cdot x_2+c=y_2^{\prime}+v_{y_2} \\&d \cdot y_2+\mathrm{f} \cdot x_2+g=x_2^{\prime}+v_{x_2}\\&a \cdot y_3+b \cdot x_3+c=y_3^{\prime}+v_{y_3} \\&d \cdot y_3+\mathrm{f} \cdot x_3+g=x_3^{\prime}+v_{x_3} \\&a \cdot y_4+b \cdot x_4+c=y_4^{\prime}+v_{y_4} \\&d \cdot y_4+\mathrm{f} \cdot x_4+g=x_4^{\prime}+v_{x_4} \end{aligned}$$Ako zapišemo u matričnom obliku bilo bi:$$\left[\begin{array}{cccccc}y_1 & x_1 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & y_1 & x_1 & 1 \\y_2 & x_2 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & y_2 & x_2 & 1 \\y_3 & x_3 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & y_3 & x_3 & 1 \\y_4 & x_4 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & y_4 & x_4 & 1\end{array}\right] \cdot\left[\begin{array}{c}a\\b \\c \\d \\f \\g\end{array}\right]=\left[\begin{array}{c}y'_1 \\x'_1 \\y'_2 \\x'_2 \\y'_3 \\x'_3 \\y'_4 \\x'_4\end{array}\right]+\left[\begin{array}{c}v_{y_1} \\v_{x_1} \\v_{y_2} \\v_{x_2} \\v_{y_3} \\v_{x_3} \\v_{y_4} \\v_{x_4}\end{array}\right]$$U matričnom obliku možemo zapisati i na kraći način:$$\mathbf{A} t=L+v$$Ako tražimo rešenje da na svim tačkama odstupanje bude minimalno, imaćemo:$$\begin{aligned}& v^{T} v = minimum \\& \mathbf{A}^{T}(L-\mathbf{A} \hat{t})=\mathbf{0} \\& \mathbf{A}^{T} L=\left(\mathbf{A}^{T} \mathbf{A}\right) \hat{t} \\& \hat{t}=\left(\mathbf{A}^{T} \mathbf{A}\right)^{-1}\left(\mathbf{A}^{T} L\right)\end{aligned}$$Uzmimo primer iz rada „Osvrt na afinu transformaciju“ [@lapaine1994osvrt] i sračunajmo transformacione parametre za zadate identične tačke.| **y** | **x** | **y'** | **x'** ||:--------------:|:--------:|:--------------:|:--------------:|| -7117.74 | 19198.04 | 398349.19 | 934806.24 || -2749.10 | 18090.87 | 402701.38 | 933635.66 || -3778.81 | 15540.10 | 401634.66 | 931100.14 || -8051.19 | 15535.25 | 397362.55 | 931157.39 |```{r aff, message=FALSE, warning=FALSE}options(digits =12)y =c(-7117.74 , -2749.1, -3778.81, -8051.19)x =c(19198.04, 18090.87, 15540.1, 15535.25)y1 =c( 398349.19, 402701.38, 401634.66, 397362.55 )x1 =c( 934806.24, 933635.66, 931100.14, 931157.39 )A =rbind(c( y[1] , x[1], 1, 0,0,0 ) ,c(0,0,0, y[1] , x[1], 1 ),c( y[2] , x[2], 1, 0,0,0 ),c(0,0,0, y[2] , x[2], 1),c( y[3] , x[3], 1, 0,0,0 ),c(0,0,0, y[3] , x[3], 1 ),c( y[4] , x[4], 1, 0,0,0 ),c(0,0,0, y[4] , x[4], 1 ))L =rbind(y1[1] , x1[1], y1[2] , x1[2], y1[3] , x1[3], y1[4] , x1[4])tkoef =round(solve(t(A)%*%A)%*%t(A)%*%L, digits=8)# transformacioni parametritkoef# Transformacija koordinata yxr = A%*%tkoefyr = yxr[seq(1,7,by=2)]xr = yxr[seq(2,8,by=2)]print(round( cbind(yr,xr), 2) )# suma kvadrata odstupanja sum ( (L-yxr)^2 )```U geodeziji je uobičajeno da afinu transformaciju definišemo na sledeći način kad imamo definisan jedan parametar razmere (skaliranja) koji je isti za obe koordinatne ose, a koordinate su zarotirane za ugao $\theta$:$$\left[\begin{array}{l}E^{\prime} \\N^{\prime}\end{array}\right]=\mu \left[\begin{array}{ll}\cos \theta & -\sin \theta \\\sin \theta & \cos \theta \\\end{array}\right]\left[\begin{array}{l}E\\N\end{array}\right]+\left[\begin{array}{l}t_E \\t_N\end{array}\right]$$Gde je $\mu$ parametar razmere (skaliranja), $\theta$ ugao rotacije i $t_E, t_N$ parametri translacije.Slično bismo dalje mogli prikazati afinu transformaciju u 3D prostoru:$$\left[\begin{array}{l}x^{\prime} \\y^{\prime} \\z^{\prime}\end{array}\right]=\left[\begin{array}{lll}a_{11} & a_{12} & a_{13} \\a_{21} & a_{22} & a_{23} \\a_{31} & a_{32} & a_{33}\end{array}\right]\left[\begin{array}{l}x \\y \\z\end{array}\right]+\left[\begin{array}{l}b_1 \\b_2 \\b_3\end{array}\right], \operatorname{det}\left[a_{i j} \right] \neq 0$$Gde bi se mogla izdvojiti matrica transformacije:$$A_b=\left[\begin{array}{cccc}a_{11} & a_{12} & a_{13} & b_1 \\a_{21} & a_{22} & a_{23} & b_2 \\a_{31} & a_{32} & a_{33} & b_3 \\0 & 0 & 0 & 1\end{array}\right]$$U `proj` biblioteke implementiran je veliki broj transformacija. Pogledati detalje na url adresi: <https://proj.org/en/9.2/operations/transformations/affine.html>):$$\left[\begin{array}{c}X \\Y \\Z \\T\end{array}\right]^{\text {dest }}=\left[\begin{array}{c}\text { xoff } \\\text { yoff } \\\text { zoff } \\\text { toff }\end{array}\right]+\left[\begin{array}{cccc}s 11 & s 12 & s 13 & 0 \\s 21 & s 22 & s 23 & 0 \\s 31 & s 32 & s 33 & 0 \\0 & 0 & 0 & \text { tscale }\end{array}\right]\left[\begin{array}{c}X \\Y \\Z \\T\end{array}\right]^{\text {source }}$$Iz dokumentacije `sf::vignette("sf3")`[@sf] primer ilustruje afinu transformaciju. Sve države Severne Karoline su zarotirane za $90^o$ oko centroida svakog poligona i skalirane faktorom 0.75, @fig-aff.```{r fig-aff, message=FALSE, warning=FALSE}#| fig.cap: "Afina transformacija, primer iz sf paketa, vinjeta sf3."#| out.width: 100%#| code-fold: falselibrary(sf)# matrica rotacijerot =function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)nc =st_read(system.file("shape/nc.shp", package="sf"), quiet =TRUE)ncg =st_geometry(nc)plot(ncg, border ='grey')cntrd =st_centroid(ncg)ncg2 = (ncg - cntrd) *rot(pi/2) * .75+ cntrdplot(ncg2, add =TRUE)plot(cntrd, col ='red', add =TRUE, cex = .5)```U okviru `proj` bibliteke implementiran je veliki broj transformacija. Pogledati detalje na url adresi:<https://proj.org/en/9.2/operations/transformations/index.html>.
