Larus audouinii - obrada GPS podataka

Postavke

# minimalni broj lokacija, foragingtrips() funkcija će filtrirati kraće tripove:
opt_minp <- 3
# radijus oko gnijezda za definiranje "kolonije" (u metrima):
opt_radius <- 500
# definiranje CRS-ova
wgs84 <- crs("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
htrs96 <- crs("+proj=tmerc +lat_0=0 +lon_0=16.5 +k=0.9999 \
              +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
# definiranje mjesta gniježđenja (središte kolonije)
nest <- SpatialPoints(list(longitude = 549608, latitude = 4735945), proj4string = htrs96)

# definiranje spola
sextab <- data.frame(GpsID = c("CROG01", "CROG02", "CROG03", "CROG04", "CROG05"), 
                     Sex = as.factor(c("Female", "Female", "Male", "Male", "Female")))

Unos podataka

Granice država

# dohvaćanje podloga
# funkcija raster::getData() dohvaća s GADM online baze državne granice:
hr <- getData(name = "GADM", country = "HRV", level = 0, path = 'data/')
it <- getData(name = "GADM", country = "ITA", level = 0, path = 'data/')
al <- getData(name = "GADM", country = "ALB", level = 0, path = 'data/')
mn <- getData(name = "GADM", country = "BIH", level = 0, path = 'data/')
bh <- getData(name = "GADM", country = "MNE", level = 0, path = 'data/')
sl <- getData(name = "GADM", country = "SVN", level = 0, path = 'data/')

# transformacija CRS-a za HR granicu
hr_proj <- spTransform(hr, CRSobj = htrs96)

# spajanje pojedinačnih država u jedan sloj
drzave <- rbind(hr, it, al, mn, bh, sl) %>% spTransform(CRSobj = htrs96)
plot(drzave)

# manje precizne granice za cijeli svijet
svijet <- readOGR(dsn = 'data/worldborders/', layer = 'TM_WORLD_BORDERS-0.3')

## OGR data source with driver: ESRI Shapefile 
## Source: "data/worldborders/", layer: "TM_WORLD_BORDERS-0.3"
## with 246 features
## It has 11 fields

Ostale podloge za pribaviti

• Batimetrija
• Temperatura mora
• Prim. produkcija (klorofil a)
• Kretanje ribarskih brodica

GPS točke

Dohvat podataka

GPS podatke s Ecotone servera direktno dohvaća funkcija get_ecotone_data() koja se nalazi u odvojenoj skripti s funkcijama telemetry_f.R. Zatim odvojene mjesece posprema u jednu tablicu, tab, koju funkcijom mkspat() pretvara u R-ov prostorni format (sptab).

# username i šifra za pristup ecotone website-u
auth <- readLines("data/auth_laraud.txt")

# funkcija get_ecotone_data() dohvaća sa ecotone servera posljednju verziju podataka za određeni mjesec i godinu:
suppressMessages({
mj05 <- get_ecotone_data(2017, 5, "crogull", auth[1], auth[2])
mj06 <- get_ecotone_data(2017, 6, "crogull", auth[1], auth[2])
mj07 <- get_ecotone_data(2017, 7, "crogull", auth[1], auth[2])
mj08 <- get_ecotone_data(2017, 8, "crogull", auth[1], auth[2])
})

tab <- bind_rows(mj05, mj06, mj07, mj08) %>% 
  rename(Longitude = Longtitude) # tipfeler u Ecotone formatu
tab$GpsID <- as.factor(tab$GpsDescription)

# mkspat() funkcija od obične tablice radi prostornu tablicu, tj. spatialpoints dataframe
# movestats() fja dodaje point-to-point udaljenosti, vremensku razliku i brzinu
sptab <- mkspat(tab, crs = htrs96) %>% 
  movestats() %>% 
  left_join(sextab, by = "GpsID") # dodaje spol, prethodno zadan u sextab varijabli
plot(mkspat(tab, crs = wgs84), col = sptab$GpsID, pch = "o")
plot(svijet, add = TRUE)

Obrada podataka

Kategorizacija prema prisutnosti na gnjezdištu

Dodan je logički vektor sptab$AtColony koji govori je li pojedina GPS točka ‘na koloniji’ (TRUE) ili nije (FALSE). Radijus kolonije definiran je u varijabli opt_radius (vidi “Postavke”).

# dodavanje logičkog vektora za prisutnost na koloniji
sptab$AtColony <- spatfilter(input = sptab, excl_geom = nest, radius = opt_radius)
{
  plot.c(drzave, radius = 1, col = "grey")
  plot(gBuffer(nest, width = opt_radius), add = TRUE)
  plot(sptab[sptab$AtColony,], add = TRUE)
  plot(sptab[!sptab$AtColony,], add = TRUE, col = rgb(1, 0, 0))
}

Filtriranje prema periodu gniježđenja

U varijabli period su zadani početak (tj. vrijeme označavanja) i kraj gniježđenja.

# definiranje perioda gniježđenja za svaku jedinku
period <- tibble(GpsID = as.factor(c("CROG01", "CROG02", "CROG03", "CROG04", "CROG05")),
                 # ovdje je definiran početak sezone, tj. za 2017. trenutak postavljanja trackera 
                 pocetak = as.POSIXct(c("2017-05-02 11:00:00",
                                        "2017-05-02 11:00:00",
                                        "2017-05-03 13:00:00",
                                        "2017-05-03 13:00:00",
                                        "2017-05-03 13:00:00"), 
                                      format = "%Y-%m-%d %H:%M:%S", 
                                      tz = "UTC"),
                 # ovdje je definiran kraj gniježđenja, tj. trenutak od kojeg se iz ponašanja ptice
                 # može zaključiti da više ne gnijezdi ("ručno" utvrđeno uvidom u podatke):
                 kraj = as.POSIXct(c("2017-06-04 03:00:00",
                                     "2017-07-09 12:00:00",
                                     "2017-06-30 10:00:00",
                                     "2017-06-08 07:30:00",
                                     "2017-06-25 07:00:00"), 
                                   format = "%Y-%m-%d %H:%M:%S", 
                                   tz = "UTC"))

# filtriranje prema periodu gniježđenja
sptab_g <- sptab %>%
  left_join(period, by = "GpsID") %>% 
  filter(GPSTime > pocetak & GPSTime < kraj)

{
  plot.p(drzave)
  plot(sptab_g, col = sptab$GpsID, pch = '*', add = TRUE)
}

Prepoznavanje izlazaka s kolonije

Izlasci s kolonije, ili “foraging trips”, su definirani kao skup uzastopnih GPS lokacija koje imaju preko neke minimalne vrijednosti zračne udaljenosti (opt_radius: 500) od gnijezda (varijabla nest). Funkcija foragingtrips() pronalazi takve uzastopne GPS točke i numerira ih prema izlasku u stupac TripID.

# funkcija za numeriranje odlazaka s kolonije
# ulaz: dataframe sa $AtColony stupcem (i.e. output spatfilter() funkcije)
foragingtrips <- function(sp_table, minpoints = 3){
  sp_table <- sp_table %>% 
    mutate(TripID = 1 + c(0, cumsum(abs(diff(sp_table$AtColony))))) %>% 
    filter(!AtColony) %>% 
    mutate(TripID = as.numeric(as.factor(TripID)))
  
  trippoints <- sp_table %>% 
    as.tibble() %>% 
    group_by(TripID) %>% 
    summarise(TripPoints = n())

  sp_table <- sp_table %>% 
    left_join(trippoints, by = "TripID") %>% 
    filter(TripPoints >= minpoints)
  # izlaz: dataframe sa TripID i TripPoints stupcima
  # TripID: redni broj za "foraging trip"
  # TripPoints: ukupni broj GPS lokacija za taj foraging trip (za potrebe filtriranja)
  return(sp_table)
}

Stupac TripPoints za svaku točku daje broj ukupnih GPS točaka zabilježenih u tom izlasku.

sptab_g %>%
  filter(GpsID == "CROG02") %>% 
  foragingtrips(., minpoints = opt_minp) %>% 
  dplyr::select(GpsID, TripID, TripPoints) %>%
  as.data.frame %>% 
  head(., 25)

##     GpsID TripID TripPoints coords.x1 coords.x2
## 1  CROG02      1          5  538145.6   4753948
## 2  CROG02      1          5  537351.1   4753425
## 3  CROG02      1          5  538635.5   4754186
## 4  CROG02      1          5  538939.0   4754173
## 5  CROG02      1          5  549333.5   4744703
## 6  CROG02      2         10  559111.5   4749205
## 7  CROG02      2         10  575384.0   4752147
## 8  CROG02      2         10  583729.1   4756628
## 9  CROG02      2         10  590444.8   4751069
## 10 CROG02      2         10  590441.3   4752056
## 11 CROG02      2         10  590368.0   4751322
## 12 CROG02      2         10  590040.4   4751594
## 13 CROG02      2         10  576264.4   4750320
## 14 CROG02      2         10  568793.2   4740445
## 15 CROG02      2         10  555025.3   4736239
## 16 CROG02      3         27  553916.2   4742001
## 17 CROG02      3         27  558637.6   4754166
## 18 CROG02      3         27  563101.6   4756570
## 19 CROG02      3         27  570071.4   4753359
## 20 CROG02      3         27  575090.8   4752814
## 21 CROG02      3         27  579731.0   4754177
## 22 CROG02      3         27  579710.4   4754079
## 23 CROG02      3         27  590056.6   4751187
## 24 CROG02      3         27  580274.0   4754295
## 25 CROG02      3         27  568629.7   4761426

# numeriranje odlazaka s kolonije
g1 <- foragingtrips(filter(sptab_g, GpsID == "CROG01"), minpoints = opt_minp)
g2 <- foragingtrips(filter(sptab_g, GpsID == "CROG02"), minpoints = opt_minp)
g3 <- foragingtrips(filter(sptab_g, GpsID == "CROG03"), minpoints = opt_minp)
g4 <- foragingtrips(filter(sptab_g, GpsID == "CROG04"), minpoints = opt_minp)
g5 <- foragingtrips(filter(sptab_g, GpsID == "CROG05"), minpoints = opt_minp)

# izbacivanje prvog (prekratkog) tripa za g5
g5 <- g5 %>% filter(TripID != 1)

# generiranje trip objekata za sve track-ove
g1_trip <- trip(g1, TORnames = c("GPSTime", "TripID"))
g2_trip <- trip(g2, TORnames = c("GPSTime", "TripID"))
g3_trip <- trip(g3, TORnames = c("GPSTime", "TripID"))
g4_trip <- trip(g4, TORnames = c("GPSTime", "TripID"))
g5_trip <- trip(g5, TORnames = c("GPSTime", "TripID"))

# kreiranje "trajectory" objekta
g1_traj <- as.ltraj(xy = coordinates(g1), date = g1$GPSTime, id = g1$TripID)
g2_traj <- as.ltraj(xy = coordinates(g2), date = g2$GPSTime, id = g2$TripID)
g3_traj <- as.ltraj(xy = coordinates(g3), date = g3$GPSTime, id = g3$TripID)
g4_traj <- as.ltraj(xy = coordinates(g4), date = g4$GPSTime, id = g4$TripID)
g5_traj <- as.ltraj(xy = coordinates(g5), date = g5$GPSTime, id = g5$TripID)

Opis kretanja

Podaci su grupirani prema pticama i izlascima s kolonije (foraging trips, u daljnjem tekstu “izlasci”) te su izračunati sljedeći podaci:

• pocetak - vrijeme odlaska s kolonije.
• kraj - vrijeme povratka na koloniju.
• pocetak_mj - mjesec u kojem je jedinka otišla s kolonije, kao nominalna varijabla.
• kraj_mj - mjesec u kojem se jedinka vratila na koloniju, kao nominalna varijabla.
• l_p2p_total - ukupna prijeđena udaljenost tijekom izlaska \([km]\).
• l_nest_mean - srednja zračna udaljenost od gnijezda tijekom izlaska \([km]\).
• l_nest_max - najveća zračna udaljenost od gnijezda tijekom izlaska \([km]\).
• nest_max_x, nest_max_y - x i y koordinate najudaljenije točke od gnijezda tijekom izlaska.
• srednja_brzina - srednja brzina kretanja tijekom izlaska \([m/s]\).

Cijela tablica je snimljena u laraud_trips.csv.

# ponovno spaja podatke za pojedine galebove u jednu tablicu
g <- rbind(g1, g2, g3, g4, g5)

gt <- g %>% 
  as.tibble() %>% 
  group_by(GpsID, TripID) %>% 
  summarise(pocetak = min(GPSTime), 
            kraj = max(GPSTime), 
            pocetak_mj = month(min(GPSTime), label = TRUE), 
            kraj_mj = month(max(GPSTime), label = TRUE),
            l_p2p_total = sum(Distance) / 1000, 
            l_nest_mean = mean(Dist_from_nest) / 1000,
            l_nest_median = median(Dist_from_nest) / 1000,
            l_nest_max = max(Dist_from_nest) / 1000,
            nest_max_x = LonProj[[which.max(Dist_from_nest)]],
            nest_max_y = LatProj[[which.max(Dist_from_nest)]],
            # trajanje = sum(Timediff) / 3600, # stara metoda za izračun
            srednja_brzina = mean(Speed, na.rm = TRUE)) %>% 
  mutate(trajanje = difftime(kraj, pocetak, units = 'hours')) %>% 
  left_join(sextab, by = "GpsID")

Ovom metodom su u izlascima (“foraging trips”) identificirane 7938 GPS točke, od ukupno zabilježenih 16197 GPS točaka.

Isti podaci su izračunati i skupno po jedinci i prikazani u donjoj tablici, iz koje je vidljivo:

1. Po varijabli l_nest_mean vidi se da postoji jasna razlika između ponašanja jedinki CROG02 i CROG04 te svih ostalih s obzirom na zadržavanje u blizini kolonije.
2. Po varijabli l_p2p_total vidi se da ne postoji velika razlika između ukupne udaljenosti koju su prešle ptice koje su se tijekom vremena istraživanja zadržavale u blizini kolonije (CROG02, CROG04) i ostalih ptica.

gi <- g %>% 
  as.tibble() %>% 
  group_by(GpsID) %>% 
  summarise(pocetak = min(GPSTime), 
            kraj = max(GPSTime), 
            pocetak_mj = month(min(GPSTime), label = TRUE), 
            kraj_mj = month(max(GPSTime), label = TRUE),
            l_p2p_total = sum(Distance) / 1000, 
            l_nest_mean = mean(Dist_from_nest) / 1000,
            l_nest_median = median(Dist_from_nest) / 1000,
            l_nest_max = max(Dist_from_nest) / 1000,
            nest_max_x = LonProj[[which.max(Dist_from_nest)]],
            nest_max_y = LatProj[[which.max(Dist_from_nest)]],
            # trajanje = sum(Timediff) / 3600, # stara metoda za izračun
            srednja_brzina = mean(Speed, na.rm = TRUE)) %>% 
  mutate(trajanje = difftime(kraj, pocetak, units = 'hours')) %>% 
  left_join(sextab, by = "GpsID")

gi %>% 
  dplyr::select(-pocetak_mj, -kraj_mj, -nest_max_x, -nest_max_y, -Sex) %>% 
  knitr::kable()

GpsID	pocetak	kraj	l_p2p_total	l_nest_mean	l_nest_median	l_nest_max	srednja_brzina	trajanje
CROG01	2017-05-02 14:00:11	2017-06-03 12:30:24	3572.563	186.45974	200.92648	220.55897	1.725691	766.5036 hours
CROG02	2017-05-02 11:00:14	2017-07-09 10:00:19	7283.771	29.43883	32.90945	50.28673	1.590056	1631.0014 hours
CROG03	2017-05-04 04:00:13	2017-06-29 20:30:26	6853.602	145.65600	213.51673	246.95580	1.994603	1360.5036 hours
CROG04	2017-05-03 14:00:15	2017-06-08 07:00:15	4690.498	48.54801	25.82479	219.89667	2.376933	857.0000 hours
CROG05	2017-05-06 07:00:21	2017-06-25 06:30:13	6346.502	161.84477	219.55259	255.92294	2.444955	1199.4978 hours

Pretvaranje u kretnice

Korišteni paketi: adehabitatLT i trip.

# generiranje trip objekata za sve track-ove
g1_trip <- trip(g1, TORnames = c("GPSTime", "TripID"))
g2_trip <- trip(g2, TORnames = c("GPSTime", "TripID"))
g3_trip <- trip(g3, TORnames = c("GPSTime", "TripID"))
g4_trip <- trip(g4, TORnames = c("GPSTime", "TripID"))
g5_trip <- trip(g5, TORnames = c("GPSTime", "TripID"))

# kreiranje "trajectory" objekta
g1_traj <- as.ltraj(xy = coordinates(g1), date = g1$GPSTime, id = g1$TripID)
g2_traj <- as.ltraj(xy = coordinates(g2), date = g2$GPSTime, id = g2$TripID)
g3_traj <- as.ltraj(xy = coordinates(g3), date = g3$GPSTime, id = g3$TripID)
g4_traj <- as.ltraj(xy = coordinates(g4), date = g4$GPSTime, id = g4$TripID)
g5_traj <- as.ltraj(xy = coordinates(g5), date = g5$GPSTime, id = g5$TripID)
{
  plot(g1_trip, main = "CROG01")
  lines(drzave)
  lines(g1_trip)
  
  plot(g2_trip, main = "CROG02")
  lines(drzave)
  lines(g2_trip)
  
  plot(g3_trip, main = "CROG03")
  lines(drzave)
  lines(g3_trip)
  
  plot(g4_trip, main = "CROG04")
  lines(drzave)
  lines(g4_trip)
  
  plot(g5_trip, main = "CROG05")
  lines(drzave)
  lines(g5_trip)
}  

Opis kretanja prema jedinki i spolu

Korišten paket ggplot2.

p <- ggplot(gt) + 
  labs(x = "Jedinka")
  
p + geom_boxplot(aes(GpsID, l_p2p_total)) + labs(y = "Ukupna duljina puta tijekom izlaska [m]")

p + geom_boxplot(aes(GpsID, srednja_brzina)) + labs(y = "Srednja brzina tijekom izlaska [m/s]")

p + geom_boxplot(aes(GpsID, l_nest_mean)) + labs(y = "Srednja zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(GpsID, l_nest_median)) + labs(y = "Srednja zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(GpsID, l_nest_max)) + labs(y = "Maksimalna zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(GpsID, as.numeric(trajanje))) + labs(y = "Ukupno trajanje izlaska [h]")

p <- ggplot(gt) +
  labs(x = "Spol")

p + geom_boxplot(aes(Sex, l_p2p_total)) + labs(y = "Ukupna duljina puta tijekom izlaska [m]")

p + geom_boxplot(aes(Sex, srednja_brzina)) + labs(y = "Srednja brzina tijekom izlaska [m/s]")

p + geom_boxplot(aes(Sex, l_nest_mean)) + labs(y = "Srednja zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(Sex, l_nest_median)) + labs(y = "Srednja zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(Sex, l_nest_max)) + labs(y = "Maksimalna zračna udaljenost od gn. tijekom izlaska [m]")

p + geom_boxplot(aes(Sex, as.numeric(trajanje))) + labs(y = "Ukupno trajanje izlaska [h]")

Vrijeme odlaska i povratka na koloniju, prema vremenskim kategorijama:

• p_kat - vrijeme odlaska
• k_kat - vrijeme povratka

gt$p_kat <- NA
gt$k_kat <- NA

# timecat() kategorizira datetime varijablu u 4 kategorije
gt$p_kat <- timecat(gt$pocetak, 
                    cutpoints = c(5, 11, 17, 23), # u UTC vremenskoj zoni
                    categories = c("01-07 h",
                                   "07-13 h",
                                   "13-19 h",
                                   "19-01 h"))
gt$k_kat <- timecat(gt$kraj, 
                    cutpoints = c(5, 11, 17, 23),
                    categories = c("01-07 h",
                                   "07-13 h",
                                   "13-19 h",
                                   "19-01 h"))

gggt <- gt %>% 
  dplyr::select(GpsID, TripID, pocetak, kraj, Odlazak = p_kat, Povratak = k_kat, Sex) %>% 
  gather(oznaka, vrijeme, Odlazak:Povratak)

p <- ggplot(data = gggt, aes(x = vrijeme, fill = oznaka))
p + geom_bar(stat = "count", position = "dodge") + facet_grid(~Sex) + 
  labs(y = "frekvencija")

p <- ggplot(data = gggt, aes(x = vrijeme, fill = Sex))
p + geom_bar(stat = "count", position = "dodge") + facet_grid(~oznaka) + 
  labs(y = "frekvencija")

Aktivnost danju/noću

Dan i noć definirani su zavisno od izlaska i zalaska sunca prema datumu i lokaciji koristeći NOAA tablicu efemerida (funkcija sunriset iz paketa maptools).

dayornight <- function(latitude, longitude, datetime) {
  library(maptools)
  location <- matrix(c(longitude, latitude), nrow=1)
  sunrise <- sunriset(location, datetime, direction = "sunrise", POSIXct.out = TRUE)
  sunset <- sunriset(location, datetime, direction = "sunset", POSIXct.out = TRUE)
  # print(sunrise)
  # print(sunset)
  if (datetime < sunrise$time | datetime > sunset$time) {return("night")}
  if (datetime >= sunrise$time & datetime <= sunset$time) {return("day")}
}

Lokacija na kojoj je jedinka aktivna definirana je kao GPS točka minimalno 300 m udaljena od prethodne zabilježene GPS točke, što eliminira GPS ‘drift’ i vrlo vjerojatno upućuje na aktivno kretanje. Potencijalni problem za ovaj način definiranja može biti odmaranje na moru u uvjetima s prisutnom površinskom morskom strujom.

tab2 <- as_tibble(sptab)
tab2$daynight <- NA

for (i in 1:nrow(tab2)) {
  tab2[i,]$daynight <- dayornight(latitude = tab2[i,]$Latitude,
                                  longitude = tab2[i,]$Longitude,
                                  datetime = tab2[i,]$GPSTime)
}

## Checking rgeos availability: TRUE

tab2$isactive <- FALSE
tab2[tab2$Distance >= 300,]$isactive <- TRUE

p4 <- ggplot(data = tab2, aes(x = isactive, y = (..count..)/sum(..count..)))
p4 + geom_bar(stat = "count", position = "dodge") + facet_grid(~daynight) +
  labs(y = "rel. frequency", x = "activity")

Iz gornjeg grafa je vidljivo da je razina aktivnosti znatno veća tijekom dana, što potvrđuje i hi-kvadrat test:

conttab <- tab2 %>% dplyr::select(daynight, isactive) %>% table
hikv <- chisq.test(conttab)

hikv

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  conttab
## X-squared = 1360.3, df = 1, p-value < 2.2e-16

hikv$observed

##         isactive
## daynight FALSE TRUE
##    day    3866 6133
##    night  4245 1953

hikv$expected

##         isactive
## daynight    FALSE     TRUE
##    day   5007.217 4991.783
##    night 3103.783 3094.217

Za vremensku analizu potrebno je isključiti točke za koje je Timediff varijabla jako dugačka, jer ako je prošlo više sati ili dana, ne možemo znati je li ptica bila aktivna za vrijeme dana ili noći. Ovdje je 3600 sekundi uzeto kao cutoff vrijednost.

activehours <- tab2 %>% 
  filter(Timediff < 3600) %>% 
  group_by(isactive, daynight) %>% 
  summarise(activehours = sum(Timediff)/3600)

sati_noc <- activehours %>% filter(daynight == "night") %>% .$activehours %>% sum
akt_sati_noc <- activehours %>% filter(daynight == "night", isactive) %>% .$activehours %>% sum
akt_sati_noc/sati_noc * 100

## [1] 28.08149

sati_dan <- activehours %>% filter(daynight == "day") %>% .$activehours %>% sum
akt_sati_dan <- activehours %>% filter(daynight == "day", isactive) %>% .$activehours %>% sum
akt_sati_dan/sati_dan * 100

## [1] 59.46125

Dakle noću su ptice bile aktivne 28.0814883% vremena, a danju 59.4612461% vremena.

hikv2 <- activehours %>% 
  xtabs(activehours ~ daynight + isactive, .) %>% 
  chisq.test()

hikv2

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  .
## X-squared = 660.31, df = 1, p-value < 2.2e-16

hikv2$observed

##         isactive
## daynight     FALSE      TRUE
##    day   1755.0281 2574.2319
##    night 1964.4861  767.0583

hikv2$expected

##         isactive
## daynight    FALSE     TRUE
##    day   2280.582 2048.678
##    night 1438.932 1292.612

ggplot(data = activehours, aes(x = isactive, y = activehours)) + 
  geom_bar(stat = "identity") + 
  facet_grid(~daynight)

Kretanje ribarskih brodica

Učitani su podaci o kretanju ribarskih brodica za period gniježđenja galebova.

brod <- read.csv(file = "data/brodovi/brodovi.csv", sep = ";", header = TRUE) %>%
  tbl_df() %>% 
  mutate(time = dmy_hm(GPS_DTIME)) %>% 
  rename(Latitude = WGS84_LAT, Longitude = WGS84_LONG) %>% 
  nodups %>% 
  mkspat(., crs = htrs96) %>% 
  mutate(DateTime = dmy_hm(GPS_DTIME))
plot(brod, pch = "."); lines(hr)

Područje za ovu analizu definirano je kao minimalni konveksni poligon koji obuhvaća 99% GPS lokacija ribarskih brodova iz dostupnog seta podataka, iz kojeg je izrezano kopno.

brodmcp <- brod %>% mutate(id = 1) %>% dplyr::select(id) %>% mcp(percent = 100)
razlika <- gDifference(brodmcp, drzave)
plot(brodmcp)
lines(drzave, col = rgb(0.4, 0.4, 0.4))
lines(razlika, col = rgb(1, 0, 0))

Unutar područja analize nasumično je generirano 5000 točaka (kontrola).

randomsample <- spsample(x = razlika, n = 5000, type = "random")
plot(randomsample)
plot(drzave, add = TRUE)

Za svaku GPS lokaciju galeba izračunate su sljedeće mjere: * dist_to_ship - minimalna udaljenost od lokacija brodica u rasponu ±7200 sekundi od trenutka bilježenja GPS točke * dist_to_random - minimalna udaljenost od kontrole, tj. 5000 točaka nasumično raspoređenih unutar područja analize

Budući da sad ne gledamo gledamo prosječnu nego minimalnu udaljenost, ovaj oblik kontrole s 5000 nasumičnih točaka po meni nema smisla, jer veličina uzorka utječe na vrijednost – što je veći n za kontrolu, veća je šansa da će se naći točka blizu GPS lokacije galeba.

brodovi_sad <- function(brod, datetime, window) {
  diff_temp <- abs(difftime(brod$DateTime, datetime, units = "s"))
  return(brod[diff_temp < window,])
}

p2p_dist <- function(x1, y1, x2, y2) {
  dx <- x1 - x2
  dy <- y1 - y2
  return(sqrt(dx^2 + dy^2))
}

min_dist <- function(x, y, comparison_table) {
  distances <- array()
  for (i in 1:nrow(comparison_table)) {
    distances[i] <- p2p_dist(x, y, comparison_table[i,1], comparison_table[i,2])
  }
  return(min(distances))
}

sptab_br <- sptab[!is.na(over(sptab, razlika)),]
sptab_br$dist_to_ship <- NA
sptab_br$dist_to_random <- NA
for (i in 1:nrow(sptab_br)) {
  comparison_table <- brodovi_sad(brod, 
                                  sptab_br$GPSTime[i], 
                                  window = 7200)@coords
  
  if(nrow(comparison_table) != 0) {
    sptab_br$dist_to_ship[i] <- min_dist(x = sptab_br$LonProj[i], 
                                      y = sptab_br$LatProj[i],
                                      comparison_table = comparison_table)
  }
  sptab_br$dist_to_random[i] <- min_dist(x = sptab_br$LonProj[i],
                                   y = sptab_br$LatProj[i],
                                   comparison_table = randomsample@coords)
}

# window = ±3600 s
sptab_br_1h <- sptab[!is.na(over(sptab, razlika)),]
sptab_br_1h$dist_to_ship <- NA
sptab_br_1h$dist_to_random <- NA
for (i in 1:nrow(sptab_br_1h)) {
  comparison_table <- brodovi_sad(brod, 
                                  sptab_br_1h$GPSTime[i], 
                                  window = 3600)@coords
  
  if(nrow(comparison_table) != 0) {
    sptab_br_1h$dist_to_ship[i] <- min_dist(x = sptab_br_1h$LonProj[i], 
                                            y = sptab_br_1h$LatProj[i],
                                            comparison_table = comparison_table)
  }
  sptab_br_1h$dist_to_random[i] <- min_dist(x = sptab_br_1h$LonProj[i],
                                            y = sptab_br_1h$LatProj[i],
                                            comparison_table = randomsample@coords)
}

# window = ±1800 s
sptab_br_30m <- sptab[!is.na(over(sptab, razlika)),]
sptab_br_30m$dist_to_ship <- NA
sptab_br_30m$dist_to_random <- NA
for (i in 1:nrow(sptab_br_30m)) {
  comparison_table <- brodovi_sad(brod, 
                                  sptab_br_30m$GPSTime[i], 
                                  window = 1800)@coords
  
  if(nrow(comparison_table) != 0) {
    sptab_br_30m$dist_to_ship[i] <- min_dist(x = sptab_br_30m$LonProj[i], 
                                             y = sptab_br_30m$LatProj[i],
                                             comparison_table = comparison_table)
  }
  sptab_br_30m$dist_to_random[i] <- min_dist(x = sptab_br_30m$LonProj[i],
                                             y = sptab_br_30m$LatProj[i],
                                             comparison_table = randomsample@coords)
}

Udaljenost od brodova

hist(sptab_br$dist_to_ship, main = "histogram udaljenosti od brodova, prozor ±7200 sekundi")

hist(sptab_br_1h$dist_to_ship, main = "histogram udaljenosti od brodova, prozor ±3600 sekundi")

hist(sptab_br_30m$dist_to_ship, main = "histogram udaljenosti od brodova, prozor ±1800 sekundi")

summary(sptab_br$dist_to_ship)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
##     13.54  18202.82  37480.01  53614.85  74359.48 231573.96      2268

summary(sptab_br_1h$dist_to_ship)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
##     34.83  19238.00  39871.22  54717.21  81250.64 231573.96      2391

summary(sptab_br_30m$dist_to_ship)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
##     34.83  19709.19  41580.58  55522.60  81754.59 231972.51      2455

Udaljenost od kontrole

hist(sptab_br$dist_to_random)

summary(sptab_br$dist_to_random)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    2.859  634.727  950.896 1035.329 1350.982 4671.624

Kernel density analiza

Iscrtane su sljedeće konture za kernel density analizu: 95%, 90%, 80%, 75%.

g2points <- SpatialPoints(sptab_g)
g2grid <- gridmaker(g2points, resolution = 500, extend = 200)
g2gridsp <- as(g2grid, "SpatialPixels")

crs(g2grid) <- htrs96
crs(g2points) <- htrs96
crs(g2gridsp) <- htrs96


kda_g2 <- kernelUD(g2points, grid = g2gridsp)

# kda_g3 <- kernelUD(g2points, grid = g2gridsp, boundary = as(hr_proj, "SpatialLines"))
# ne radi zbog ograničenja algoritma -- prerazvedena geometrija

# kda_g2a <- kernelUD(spatfilter(g2points, excl_geom = hr_proj , radius = 0), grid = g2gridsp)

plot(getverticeshr(kda_g2, percent = 95), col = rgb(1,1,1))
plot(g2points, pch = 16, col = rgb(0, 0, 0, 0.3), add = TRUE)
lines(hr_proj)
lines(getverticeshr(kda_g2, percent = 95), col = rgb(1, 0.3, 0.3))
lines(getverticeshr(kda_g2, percent = 90), col = rgb(0, 0.5, 0.5))
lines(getverticeshr(kda_g2, percent = 80), col = rgb(0.2, 0.2, 1))
lines(getverticeshr(kda_g2, percent = 75))

Kad se razmatraju sve točke iz perioda gniježđenja, rezultat nije jako informativan. Kad uzmemo samo točke iz hrvatskog teirtorijalnog mora (što bi imalo smisla ako radimo kernel analizu u svrhu određivanja zaštićenih područja) onda dobivene konture izgledaju puno smislenije.

hr_more <- readOGR(dsn = "data/worldborders/", layer = "hr_more")

## OGR data source with driver: ESRI Shapefile 
## Source: "data/worldborders/", layer: "hr_more"
## with 1 features
## It has 3 fields

crs(hr_more) <- htrs96

g2points <- SpatialPoints(sptab_g)
crs(g2points) <- htrs96

g2points <- g2points[!is.na(over(g2points, hr_more)[1]),]

g2grid <- gridmaker(g2points, resolution = 500, extend = 200)
g2gridsp <- as(g2grid, "SpatialPixels")

crs(g2grid) <- htrs96
crs(g2gridsp) <- htrs96


kda_g2 <- kernelUD(g2points, grid = g2gridsp)

# kda_g3 <- kernelUD(g2points, grid = g2gridsp, boundary = as(hr_proj, "SpatialLines"))
# ne radi zbog ograničenja algoritma -- prerazvedena geometrija

# kda_g2a <- kernelUD(spatfilter(g2points, excl_geom = hr_proj , radius = 0), grid = g2gridsp)

plot(getverticeshr(kda_g2, percent = 95), col = rgb(1,1,1))
plot(g2points, pch = 16, col = rgb(0, 0, 0, 0.3), add = TRUE)
lines(hr_proj)
lines(getverticeshr(kda_g2, percent = 95), col = rgb(1, 0.3, 0.3))
lines(getverticeshr(kda_g2, percent = 90), col = rgb(0, 0.5, 0.5))
lines(getverticeshr(kda_g2, percent = 80), col = rgb(0.2, 0.2, 1))
lines(getverticeshr(kda_g2, percent = 75))

Data manipulation & visualization

All data analysis tasks were performed in the R programming environment (R Core Team 2018). Data manipulation and visualization was performed with the packages dplyr (Wickham and Francois 2016), lubridate (Grolemund and Wickham 2011) and ggplot2 (Wickham 2009). Manipulation, analysis and visualization of spatial data was performed with the packages raster (Hijmans 2016), sp (E. J. Pebesma and Bivand 2005, Bivand, Pebesma, and Gomez-Rubio (2013)), rgeos (R. Bivand and Rundel 2017), GISTools (Brunsdon and Chen 2014), rgdal (R. Bivand, Keitt, and Rowlingson 2017) and spdplyr (Sumner 2017).

Home range estimation

For home range estimation, the kernelUD function from the adehabitatHR (Calenge 2006) package was used. The reference (ad-hoc) method was used for bandwidth estimation. Least square cross validation (LSCV) was also tried, but the algorithm did not converge and the bandwidth parameter could not be estimated with the available data.

Reference

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. http://www.asdar-book.org/.

Bivand, Roger, and Colin Rundel. 2017. Rgeos: Interface to Geometry Engine - Open Source (’Geos’). https://CRAN.R-project.org/package=rgeos.

Bivand, Roger, Tim Keitt, and Barry Rowlingson. 2017. Rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. https://CRAN.R-project.org/package=rgdal.

Brunsdon, Chris, and Hongyan Chen. 2014. GISTools: Some Further Gis Capabilities for R. https://CRAN.R-project.org/package=GISTools.

Calenge, C. 2006. “The Package Adehabitat for the R Software: Tool for the Analysis of Space and Habitat Use by Animals.” Ecological Modelling 197: 1035.

Grolemund, Garrett, and Hadley Wickham. 2011. “Dates and Times Made Easy with lubridate.” Journal of Statistical Software 40 (3): 1–25. http://www.jstatsoft.org/v40/i03/.

Hijmans, Robert J. 2016. Raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster.

Pebesma, Edzer J., and Roger S. Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.

R Core Team. 2018. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Sumner, Michael D. 2017. Spdplyr: Data Manipulation Verbs for the Spatial Classes. https://CRAN.R-project.org/package=spdplyr.

Wickham, Hadley. 2009. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. http://ggplot2.org.

Wickham, Hadley, and Romain Francois. 2016. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.