10  Rad sa rasterima i koordinatne transformacije

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.3
library(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 TRUE

bbox= ext(c(-180, 180, -90, 90) )
#(xmin, xmax, ymin, ymax)
# prazan raster 30x30 stepeni
er <- 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 piksela
er$v= 1:72
plot(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..

library(terra)
library(sf)

bbox= ext(c(-180, 180, -90, 90) )
#(xmin, xmax, ymxn, ymax)
# prazan raster 30x30 stepeni
er <- rast(bbox, resolution=30 ) 
# vrednosti piksela
er$v= 1:72
pol = 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 )

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.3

bbox= ext(c(-180, 180, -90, 90) )
#(xmin, xmax, ymin, ymax)
# prazan raster
er <- rast(bbox, resolution=30 )
# vrednosti piksela
er$v= 1:72

data(world)
granicaSRKS = st_union(world[173, 'geom' ], world[175,'geom'] )
world[173, ]$geom = granicaSRKS$geom
world = world[-175,]
# terra paket koristi svoj zapis za vektorske podatke lako je konvertovati 
# u sf objekat SpatVector
w_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 raster
er <- rast(bbox, resolution=30 )
# vrednosti piksela
er$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*1000
er_sin= project(er, "+proj=sinu", res=30*100*1000, method="near")
# zadržaćemo samo piksele koji se poklapaju sa vektorskim podacima
pol = 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 km

alt = rast("/vsicurl/https://zenodo.org/record/6211553/files/copernicus_DSM_30arcsec_mood_bbox_epsg4326.tif") 
# https://zenodo.org/record/6211553
names(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.

library(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 near
dem_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 )

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.

library(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 )

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.

library(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 cubic
dem_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 )

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.

library(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 cubicspline
dem_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 )

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.

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

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 primitiv
grf = st_point( c(20.4763281, 44.8053183) ) # lon , lat
# potom pravimo kolekciju geometrije sfc u odgovarajućoj projekciji
pts_sfc = st_sfc(grf, crs=4326)  #epsg 4326 je WGS84 long lat
# GK Srb koord. ref. sis. 3 translacije
gk3 = "+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 parametra
gk7 = "+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 UTM
pts_sfc_gk3 = st_transform(pts_sfc, crs =gk3)
pts_sfc_gk7 = st_transform(pts_sfc, crs =gk7)

# matrica koordinata iz prostornih podataka 3 translacije
st_coordinates(pts_sfc_gk3)
#                  X             Y
# [1,] 7459008.35955 4962434.20491
st_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 sekundama
plot(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.

options(digits = 12)
library(sf)
sf_proj_network(TRUE)
# [1] "https://cdn.proj.org"
b_pump <- st_read(system.file("gpkg/b_pump.gpkg", package="sf"))
# Reading layer `b_pump' from data source 
#   `C:\Users\user\AppData\Local\R\win-library\4.2\sf\gpkg\b_pump.gpkg' 
#   using driver `GPKG'
# Simple feature collection with 1 feature and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 529393.498863 ymin: 181020.577869 xmax: 529393.498863 ymax: 181020.577869
# Projected CRS: Transverse_Mercator

sf_proj_network()
# [1] TRUE
# 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=crs

bp <- sf_proj_pipelines("EPSG:27700", "EPSG:32630")
bp
# Candidate coordinate operations found:  9 
# Strict containment:     FALSE 
# Axis order auth compl:  FALSE 
# Source:  EPSG:27700 
# Target:  EPSG:32630 
# Best instantiable operation has accuracy: 1 m
# Description: Inverse of British National Grid + OSGB36 to WGS 84
#              (9) + UTM zone 30N
# Definition:  +proj=pipeline +step +inv +proj=tmerc +lat_0=49
#              +lon_0=-2 +k=0.9996012717 +x_0=400000
#              +y_0=-100000 +ellps=airy +step
#              +proj=hgridshift
#              +grids=uk_os_OSTN15_NTv2_OSGBtoETRS.tif
#              +step +proj=utm +zone=30 +ellps=WGS84
# Gridna transformacija
utm <- st_transform(b_pump, "EPSG:32630")
st_coordinates(utm)
#                  X             Y
# [1,] 698672.161451 5710795.73037

# sa helmertovom transformacijom
uk_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_pump
st_crs(b_pump_7p ) = uk_7p
utm2 <- st_transform(b_pump_7p , "EPSG:32630")
st_coordinates(utm2)
#                  X             Y
# [1,] 698672.158563 5710795.73026

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 metrima
plot(vgrid)

Slika 10.14: Vertikalni datumski grid za Veliku Britaniju u metrima, undulacije geoida.

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 } \]

10.3 Afina transformacija

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.

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:

\[ \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“ (Lapaine and Frančula 1994) 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
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 parametri
tkoef
#                 [,1]
# [1,]      0.99992012
# [2,]      0.01454191
# [3,] 405187.18493138
# [4,]     -0.01453829
# [5,]      0.99989626
# [6,] 915506.70682205

# Transformacija koordinata 
yxr = A%*%tkoef

yr = yxr[seq(1,7,by=2)]
xr = yxr[seq(2,8,by=2)]

print(round( cbind(yr,xr), 2) )
#             yr        xr
# [1,] 398349.19 934806.23
# [2,] 402701.38 933635.67
# [3,] 401634.66 931100.13
# [4,] 397362.55 931157.40

# suma kvadrata odstupanja 

sum ( (L-yxr)^2 )
# [1] 0.000173597032003

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") (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.

library(sf)
# matrica rotacije
rot = 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 + cntrd
plot(ncg2, add = TRUE)
plot(cntrd, col = 'red', add = TRUE, cex = .5)

Slika 10.15: Afina transformacija, primer iz sf paketa, vinjeta sf3.

U okviru proj bibliteke implementiran je veliki broj transformacija. Pogledati detalje na url adresi: https://proj.org/en/9.2/operations/transformations/index.html.

Lapaine, Miljenko, and Nedjeljko Frančula. 1994. “Osvrt Na Afinu Transformaciju.” Geodetski List 48 (2): 159–68.
Pebesma, Edzer. 2018. Sf: Simple Features for R. https://CRAN.R-project.org/package=sf.