########################################################################################################################### ###### DATA PREPARATION TO IDENTIFY ANTHROPOGENIC FEEDING SITES OF RED KITES BASED ON GPS TRACKING DATA ################ ########################################################################################################################### # code to replicate results in: # Predicting anthropogenic food supplementation from individual tracking data # Article DOI: 10.1111/ibi.13359 # Authors: Oppel, Steffen; Beeli, Ursin; Grüebler, Martin; van Bergen, Valentijn; Kolbe, Martin; Pfeiffer, Thomas; Scherler, Patrick # for questions contact: steffen.oppel@vogelwarte.ch ##~~~~~~~~~~~~ WARNING: this code is not self-contained! ~~~~~~~~~~~~~~~~~## ##~~~~~~~~~~~~ code only provided to understand how data in Zenodo archive were generated ~~~~~~~~~~~~~~~~~## ##~~~~~~~~~~~~ due to copyright issues the raw data cannot be openly shared by the authors ~~~~~~~~~~~~~~~~~## ####### LIBRARIES REQUIRED------------------ library(tidyverse) library(rnaturalearth) library(sf) library(amt) library(suncalc) library(dplyr, warn.conflicts = FALSE) options(dplyr.summarise.inform = FALSE) library(recurse) library(readxl) library(lubridate) library(data.table); setDTthreads(percent = 65) sf_use_s2(FALSE) # deactivating spherical geometry s2 library(move2) library(leaflet) library(units) library(foreach) library(geosphere) library(keyring) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # SET UP DOWNLOAD OF TRACKING DATA FROM MOVEBANK -------- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ####### SPECIFY THE MOVEBANK ID OF THE STUDY FOR WHICH DATA ARE SHARED MYSTUDY<-c(230545451,1356790386) MYUSERNAME<-"Steffen" movebank_store_credentials(username=MYUSERNAME, key_name = getOption("move2_movebank_key_name"), force = TRUE) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # DOWNLOAD MOVEBANK DATA AND ANIMAL INFO ---------------- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ birds1<-movebank_retrieve(study_id=MYSTUDY[1], entity_type="individual") %>% dplyr::rename(individual_id=id,bird_id=local_identifier) %>% dplyr::select(individual_id,comments, bird_id,ring_id,sex,latest_date_born) birds2<-movebank_retrieve(study_id=MYSTUDY[2], entity_type="individual") %>% dplyr::rename(individual_id=id,bird_id=local_identifier) %>% dplyr::select(individual_id,comments, bird_id,ring_id,sex,latest_date_born) birds<-bind_rows(birds1,birds2) %>% mutate(bird_id=as.numeric(as.character(bird_id))) locs1<-movebank_retrieve(study_id=MYSTUDY[1], entity_type="event", sensor_type_id="gps", timestamp_start=ymd_hms(paste(year(Sys.time())-10,"-02-15 12:00:00",sep="")), timestamp_end=ymd_hms(paste(year(Sys.time()),"-05-30 12:00:00",sep="")), progress=T) locs2<-movebank_retrieve(study_id=MYSTUDY[2], entity_type="event", sensor_type_id="gps", timestamp_start=ymd_hms(paste(year(Sys.time())-10,"-02-15 12:00:00",sep="")), timestamp_end=ymd_hms(paste(year(Sys.time()),"-05-30 12:00:00",sep="")), progress=T) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # FILTER AND COMBINE DATA ---------------- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ### MANUAL filter function to remove all locations that are within a certain time window #### turned into function to contribute to move2: https://gitlab.com/bartk/move2/-/issues/60 #' @param locs tibble. Tracking data downloaded from Movebank using \code{move2::movebank_retrieve} #' @param mintimelag numeric. Time lag (in seconds) that must have elapsed between two subsequent locations of the same individual for both locations to be retained. #' @param loc_res integer. Precision of global coordinates (in EPSG:4326) to which subsequent latitudes and longitudes are rounded when assessing whether identical locations should be filtered. mt_filter_lag<-function(locs,maxtimelag=60,loc_res=3){ out<-locs %>% group_by(individual_id) %>% mutate(prev_t=dplyr::lag(timestamp), prev_id=dplyr::lag(individual_id), prev_lat=dplyr::lag(location_lat),prev_long=dplyr::lag(location_long)) %>% mutate(dt=as.numeric(difftime(timestamp,prev_t, units="sec"))) %>% mutate(dt=dplyr::if_else(prev_id==individual_id & round(prev_lat,loc_res)==round(location_lat,loc_res) & round(prev_long,loc_res)==round(location_long,loc_res),dt,maxtimelag*2)) %>% mutate(dt=dplyr::if_else(is.na(dt),maxtimelag*2,dt)) %>% filter(dt>maxtimelag) %>% ungroup() %>% dplyr::select(timestamp,location_lat,location_long,individual_id) return(out) } filterlocs<-mt_filter_lag(bind_rows(locs1,locs2),maxtimelag=90,loc_res=1) filterlocs<-mt_filter_lag(filterlocs,maxtimelag=120,loc_res=4) trackingdata<-filterlocs %>% left_join(birds, by="individual_id") %>% mutate(tag_year=as.numeric(comments)) %>% mutate(age_cy=as.integer((timestamp-latest_date_born)/365)) %>% dplyr::select(bird_id,ring_id,sex,age_cy,timestamp,location_lat,location_long) %>% rename(long_wgs=location_long,lat_wgs=location_lat) %>% dplyr::filter(!is.na(timestamp)) %>% dplyr::filter(!is.na(lat_wgs)) %>% dplyr::filter(!is.na(bird_id)) %>% dplyr::filter(!is.na(long_wgs)) %>% filter(long_wgs<15) %>% filter(long_wgs>-10) %>% filter(lat_wgs<54) %>% filter(lat_wgs>35) %>% st_as_sf(coords = c("long_wgs", "lat_wgs"), crs=4326)%>% dplyr::mutate(long_wgs = sf::st_coordinates(.)[,1], lat_wgs = sf::st_coordinates(.)[,2]) ### LOAD INDIVIDUAL LIFE HISTORIES indseasondata<-read_excel("C:/Users/sop/OneDrive - Vogelwarte/General/DATA/Individual_life_history_2015-2023.xlsx", sheet="Individual_life_history_2015-20") %>% # updated on 3 June 2024 to include birds from 2022 dplyr::select(bird_id,ring_number,tag_year,sex_compiled, age, hatch_year) %>% rename(ring_id=ring_number) %>% mutate(hatch_year=if_else(is.na(as.numeric(hatch_year)),tag_year-3,as.numeric(hatch_year))) nestdata<-fread("C:/Users/sop/OneDrive - Vogelwarte/REKI/Analysis/NestTool/REKI/data/Basic_nest_list_2015_2022.csv") #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # DATA MANIPULATION AND PREPARATION ---------------- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # keeping only the information of relevant locations in Switzerland if("long_wgs" %in% names(trackingdata)){ trackingdata<-trackingdata %>% rename(long=long_wgs, lat=lat_wgs) } trackingdata <- trackingdata %>% filter(long>5.9) %>% filter(lat>45.8) %>% filter(long<10.6) %>% filter(lat<48) %>% filter(!is.na(timestamp)) %>% filter(!is.na(long)) %>% mutate(year_id=paste(year(timestamp),bird_id, sep="_")) %>% filter(!is.na(year_id)) dim(trackingdata) # filling in gaps in age and sex and making age uniform (in years) trackingdata <- trackingdata %>% left_join(indseasondata, by=c("bird_id","ring_id")) %>% mutate(age_cy=year(timestamp)-hatch_year) %>% mutate(sex=if_else(is.na(sex),sex_compiled,sex)) dim(trackingdata) summary(trackingdata$age_cy) head(trackingdata) # converting to metric CRS prior to estimating distances track_sf <- trackingdata %>% dplyr::select(year_id,sex,age_cy,timestamp,geometry,long,lat) %>% st_transform(crs = 3035) %>% ungroup() %>% dplyr::mutate(long_eea = sf::st_coordinates(.)[,1], lat_eea = sf::st_coordinates(.)[,2]) # Creating a track to exclude nocturnal locations (later on) track_amt <- track_sf %>% mk_track( .x = long_eea, .y = lat_eea, .t = timestamp, id = year_id, sex=sex, age_cy=age_cy, crs = 3035, order_by_ts=TRUE ) %>% time_of_day(include.crepuscule = T) %>% # if F, crepuscule is considered as night arrange(id, t_) ## clean up workspace rm(trackingdata,track_sf) gc() ### CALCULATE OTHER METRICS track_amt$step_length<-amt::step_lengths(track_amt) # include step lengths track_amt$turning_angle<-abs(amt::direction_rel(track_amt,append_last=T, full_circle = FALSE,lonlat = FALSE)) # include RELATIVE turning angles - changed from absolute track_amt$turning_angle<-ifelse(is.na(track_amt$turning_angle),0,track_amt$turning_angle) ## replace NA turning angles (caused by 0 distance) with 0 track_amt$speed<-amt::speed(track_amt) # include speed track_amt$locid<-seq_along(track_amt$t_) ### ABOVE METRICS NEED TO BE SET TO NA FOR THE LAST POSITION OF EACH INDIVIDUAL - THEY ARE NOT CALCULATED FOR YEAR_ID and therefore calculate distances between different animals lastlocs<- track_amt %>% st_drop_geometry() %>% arrange(id,t_) %>% group_by(id) %>% summarise(last=max(t_), lastloc=max(locid)) track_amt$speed[track_amt$locid %in% lastlocs$lastloc]<-NA track_amt$step_length[track_amt$locid %in% lastlocs$lastloc]<-NA track_amt$turning_angle[track_amt$locid %in% lastlocs$lastloc]<-NA head(track_amt) hist(track_amt$turning_angle*(180/pi)) range(track_amt$speed, na.rm=T) track_amt %>% filter(speed<0) track_amt %>% filter(id=="2015_139") %>% filter(yday(t_) > yday(ymd("2015-10-27"))) %>% st_drop_geometry() %>% print(n=30) track_amt %>% filter(speed>50) track_amt %>% filter(locid %in% seq(21245,21258,1)) ### FILTER OUT CRAZY SPEED LOCATIONS dim(track_amt) crazyspeedlocs<- track_amt %>% filter(speed>50) %>% st_drop_geometry() for (l in crazyspeedlocs$locid) { chunk<-track_amt %>% filter(locid %in% seq(l-5,l+5,1)) %>% st_drop_geometry() %>% select(-age_cy,-sex,-tod_) %>% mutate(prev_t=dplyr::lag(t_), post_t=dplyr::lead(t_), prev_dist=dplyr::lag(step_length),post_dist=dplyr::lead(step_length)) %>% mutate(pre_dt=as.numeric(difftime(t_,prev_t, units="sec")),post_dt=as.numeric(difftime(post_t,t_, units="sec"))) %>% rowwise() %>% mutate(dt=min(pre_dt,post_dt, na.rm=T)) %>% ungroup() ## filter the location that is (1) among the top 2 step_lengths, (2) among the top 2 turning angles, (3) min of (pre_dt and post_dt) among the bottom 2 dts sel1a<-slice_max(chunk,step_length,n=2) sel1b<-slice_max(chunk,speed,n=2) sel2<-slice_max(chunk,turning_angle,n=2) sel3<-slice_min(chunk,dt,n=2) sel3b<-slice_min(chunk,post_dt,n=1) badid<-Reduce(intersect, list(unique(sel1a$locid,sel1b$locid),sel2$locid,sel3$locid)) if(length(badid)<1){badid<-Reduce(intersect, list(sel1b$locid,sel3b$locid))} if(length(badid)<1){badid<-Reduce(intersect, list(unique(sel1a$locid,sel1b$locid),sel2$locid))} if(length(badid)==1){track_amt<-track_amt %>% filter(locid!=badid) rm(sel1,sel2,sel3,chunk,badid)} } dim(track_amt) ### RECALCULATE METRICS AFTER HAVING ELIMINATED CRAZY SPEED LOCS track_amt$locid<-seq_along(track_amt$t_) track_amt$step_length<-amt::step_lengths(track_amt) # include step lengths track_amt$turning_angle<-abs(amt::direction_rel(track_amt,append_last=T, full_circle = FALSE,lonlat = FALSE)) # include RELATIVE turning angles - changed from absolute track_amt$turning_angle<-ifelse(is.na(track_amt$turning_angle),0,track_amt$turning_angle) ## replace NA turning angles (caused by 0 distance) with 0 track_amt$speed<-amt::speed(track_amt) # include speed lastlocs<- track_amt %>% st_drop_geometry() %>% arrange(id,t_) %>% group_by(id) %>% summarise(last=max(t_), lastloc=max(locid)) track_amt$speed[track_amt$locid %in% lastlocs$lastloc]<-NA track_amt$step_length[track_amt$locid %in% lastlocs$lastloc]<-NA track_amt$turning_angle[track_amt$locid %in% lastlocs$lastloc]<-NA dim(track_amt) # RECURSIONS FOR EACH LOCATION 50 m BUFFER-------------------------------------------------- # splitting track into a list with each single id grouped to an element------ # this causes a memory limit error, so try and do it in a loop track_amt <- as.data.frame(track_amt) ### RECURSIONS IN A LOOP - takes 2 hrs -------------------------------- # calculating recursions rm(trackingdata,track_sf,buildings,forest) ### clean up workspace and memory gc() track_amt$revisits <- NA track_amt$residence_time <- NA for (i in unique(track_amt$id)){ x<-track_amt %>% filter(id==i) xr<-getRecursions(x = x[1:4], radius = 50, timeunits = "hours") track_amt$revisits[track_amt$id == i] <- xr$revisits track_amt$residence_time[track_amt$id == i] <- xr$residenceTime # CALCULATING FIRST AND LAST REVISIT AND DURATION AND TEMPORAL PERSISTENCE OF REVISITS ----------------------------------------------------------------- tempout<- xr$revisitStats %>% mutate(jday.ent=yday(entranceTime),jday.ex=yday(exitTime)) %>% group_by(coordIdx) %>% summarise(first=min(entranceTime, na.rm=T), last=max(exitTime, na.rm=T), meanFreqVisit=mean(timeSinceLastVisit, na.rm=T), n_days=length(unique(c(jday.ent,jday.ex)))) %>% mutate(TimeSpan=as.numeric(difftime(last,first,unit="days"))) %>% mutate(TempEven=n_days/TimeSpan) %>% mutate(meanFreqVisit=ifelse(is.na(meanFreqVisit),0,meanFreqVisit)) %>% ## set the frequency of visit to 0 for locations never revisited mutate(TempEven=ifelse(n_days==1,1,TempEven)) %>% ## set the evenness to 1 for locations never revisited on more than a single day select(meanFreqVisit,n_days,TimeSpan,TempEven) track_amt$meanFreqVisit[track_amt$id == i] <-tempout$meanFreqVisit track_amt$n_days[track_amt$id == i] <-tempout$n_days track_amt$TimeSpan[track_amt$id == i] <-tempout$TimeSpan track_amt$TempEven[track_amt$id == i] <-tempout$TempEven rm(tempout,x,xr) } # CALCULATING MOVING AVERAGE FOR SPEED AND ANGLE ----------------------------------------------------------------- rm(track_amt_recurse,track_amt_list,trackingdata) track_amt <- track_amt %>% arrange(id, t_) %>% group_by(id) %>% mutate(mean_speed=frollmean(speed,n=3,na.rm=T, align="center")) %>% mutate(mean_angle=frollmean(abs(turning_angle),n=5,na.rm=T, align="center")) %>% mutate(mean_speed=ifelse(is.na(mean_speed),speed,mean_speed)) %>% mutate(mean_angle=ifelse(is.na(mean_angle),turning_angle,mean_angle)) head(track_amt) dim(track_amt) ### LOADING DATA OF FORESTS AND BUILDINGS ----------------------------------------------------------------- buildings <- st_read("data/Buildings/tlm_buildings_studyareaExtra_size65_buff_50m_sf_singlepoly.shp", stringsAsFactors=FALSE) %>% st_transform(crs = 3035) forest <- st_read("data/Forest/vec25_forest_buff_20m_studyarea.shp", stringsAsFactors=FALSE) %>% st_transform(crs = 3035) %>% select(AREA) ### LOADING DATA OF ALL FEEDING STATIONS----------------------------------------------------------------- FEEDERS<- fread("data/Private_Feeders/private_feeders_upd2022.csv") %>% filter(!is.na(coordX)) %>% st_as_sf(coords = c("coordX", "coordY")) st_crs(FEEDERS) <- 21781 FEEDER_buff<- FEEDERS %>% st_transform(crs = 3035) %>% st_buffer(dist=50) %>% select(ID,type_of_food,frequency) %>% rename(feeder_id=ID) ### EXPERIMENTAL FEEDERS BY SOI ### SPATIAL DUPLICATION AT MULTIPLE TIME POINTS; SO WE TAKE THE MEDIAN FOR EACH NEST NAME EXPFEEDERS_buff<- fread("data/experimental_feeding.csv") %>% dplyr::filter(!event_id %in% c(5610,4413)) %>% ## remove 2 offending stations with duplicate coordinates that cause problems when overlaying with locations mutate(long=as.numeric(lon)) %>% filter(!is.na(long)) %>% filter(!is.na(lat)) %>% group_by(nest_name) %>% summarise(start=min(Timepoint), end=max(Timepoint), lat=mean(lat,na.rm=T),long=mean(long, na.rm=T), food_placed=median(food_placed, na.rm=T)) %>% st_as_sf(coords = c("long", "lat"), crs=4326) %>% ungroup() %>% st_transform(crs = 3035) %>% st_buffer(dist=50) %>% select(nest_name,start, end,food_placed) %>% rename(feeder_id=nest_name) st_difference(EXPFEEDERS_buff) ### FEEDING PLATFORMS PLATFORMS<- fread("data/feeding_platforms_15_16.csv") %>% filter(!is.na(x)) %>% mutate(start_date=dmy(start_date), end_date=dmy(end_date)) %>% st_as_sf(coords = c("x", "y")) st_crs(PLATFORMS) <- 21781 PLATFORMS_buff<- PLATFORMS %>% st_transform(crs = 3035) %>% st_buffer(dist=50) %>% select(Name,year,n_event,start_date,end_date) %>% rename(feeder_id=Name) PLATFORMS_buff ### LOADING DATA OF ALL NESTS----------------------------------------------------------------- nests<- nestdata %>% dplyr::select(nest_name,tree_spec,latitude,longitude) %>% st_as_sf(coords = c("longitude", "latitude")) st_crs(nests) <- 4326 nests<-nests %>% st_transform(crs = 3035) nests_buff<- nests %>% st_buffer(dist=50) nests_buff #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # SPATIAL JOINING OF TRACKING AND ENVIRONMENTAL DATA ---------------- #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## create sf object for spatial operations track_sf <- track_amt %>% st_as_sf(coords = c("x_", "y_")) # st_crs(track_sf) <- 3035 ### spatial joins with forests, buildings and feeders ## this operation bizarrely ADDs duplicate rows where polygons overlap ## need to include st_difference() for all layers to prevent this: https://gis.stackexchange.com/questions/351429/sf-st-intersection-returning-duplicate-features track_sf <- track_sf %>% st_join(st_difference(forest), join = st_intersects, left = TRUE) %>% st_join(st_difference(buildings), join = st_intersects, left = TRUE) %>% st_join(st_difference(FEEDER_buff), join = st_intersects, left = TRUE) %>% st_join(st_difference(nests_buff), join = st_intersects, left = TRUE) %>% mutate(BUILD=ifelse(is.na(build_id),0,1), NEST=ifelse(is.na(nest_name),0,1), FOREST=ifelse(is.na(AREA),0,1), FEEDER=ifelse(is.na(feeder_id),"NO","YES")) %>% ### FOR PUBLIC FEEDERS SPATIAL OVERLAP IS YES OR NO mutate(FEED_ID=feeder_id) %>% select(-feeder_id) %>% st_join(st_difference(EXPFEEDERS_buff), join = st_intersects, left = TRUE) %>% mutate(FEEDER=ifelse(is.na(feeder_id),FEEDER, ifelse(t_ %within% interval(start,end+days(30)),"YES",FEEDER))) %>% ### FOR EXPERIMENTAL FEEDERS NEED TEMPORAL OVERLAP TO WITHIN A MONTH OF Timepoint - changed to start and end of interval mutate(FEED_ID=ifelse(is.na(feeder_id),FEED_ID,feeder_id)) %>% select(-feeder_id) %>% st_join(st_difference(PLATFORMS_buff), join = st_intersects, left = TRUE) %>% mutate(FEEDER=ifelse(is.na(feeder_id),FEEDER, ifelse(t_ %within% interval(start_date,end_date+days(30)),"YES",FEEDER))) %>% ### FOR FEEDING PLATFORMS NEED TEMPORAL OVERLAP TO WITHIN A MONTH OF end time mutate(FEED_ID=ifelse(is.na(feeder_id),FEED_ID,feeder_id)) %>% select(-year,-n_event,-start_date,-end_date,-start,-end,-food_placed,-nest_name,-tree_spec,-feeder_id) %>% rename(forest_size=AREA) ### calculate distance to nearest nest ## this causes memory allocation error, so need to do it in a loop, which takes 45 min rm(track_amt,EXPFEEDERS) gc() track_sf$dist_nest <- NA for (i in unique(track_sf$id)){ x<-track_sf %>% filter(id==i) x_distances<-st_distance(x,nests) track_sf$dist_nest[track_sf$id==i] <- apply(x_distances,1,min)/1000 ### distance in km } # EXPORT THE ANNOTATED DATA FOR ALL OF SUI ----------------------------------------------------------------- track_out <- track_sf %>% st_transform(crs = 4326) %>% dplyr::mutate(long = sf::st_coordinates(.)[,1], lat = sf::st_coordinates(.)[,2]) %>% st_drop_geometry() %>% rename(year_id=id) fwrite(as.data.frame(track_out),"REKI_annotated_feeding2024_CH.csv") #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # DATA PROCESSING TO PREPARE FOR MODELLING ------------------ #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ### read in Switzerland map SUI<-readRDS("Swiss_border.rds") STUDY_AREA<-readRDS("REKI_study_area.rds") plot(STUDY_AREA) ### CALCULATE DISTANCE FROM PREDICTED FEEDING SITE TO NEAREST KNOWN SITE ############# # READ IN FEEDING LOCATIONS plot_feeders<-fread("data/Private_Feeders/private_feeders_upd2022.csv") %>% filter(!is.na(coordX)) %>% st_as_sf(coords = c("coordX", "coordY"), crs=21781) %>% st_transform(crs = 4326) %>% mutate(Type="Private") %>% select(Type) ### EXPERIMENTAL FEEDING STATIONS plot_feeders2<- fread("data/experimental_feeding.csv") %>% filter(!is.na(lon)) %>% filter(!is.na(lat)) %>% st_as_sf(coords = c("lon", "lat"), crs=4326)%>% mutate(Type="Experimental") %>% select(Type) ### FEEDING PLATFORMS - are also experimental plot_feeders3<- fread("data/feeding_platforms_15_16.csv") %>% filter(!is.na(x)) %>% mutate(start_date=dmy(start_date), end_date=dmy(end_date)) %>% st_as_sf(coords = c("x", "y"), crs=21781) %>% st_transform(crs = 4326) %>% mutate(Type="Experimental") %>% select(Type) plot_feeders<-rbind(plot_feeders,plot_feeders2,plot_feeders3) ################ FILTER DATA TO ELIMINATE FOREST AND FIELD LOCATIONS ---------------------------- ## anthropogenic feeding is unlikely in those locations track_nofor_day_build<-track_sf %>% filter(FOREST==0) %>% filter(tod_=="day") %>% filter(!(BUILD==0 & FEEDER =="NO")) %>% st_as_sf(coords = c("long", "lat"), crs = 4326) %>% st_transform(2056) %>% st_intersection(.,SUI) %>% filter(!is.na(country)) %>% ### remove all data outside of Switzerland st_transform(4326) %>% st_intersection(.,STUDY_AREA) %>% filter(!is.na(Name)) %>% ### remove all data outside of study area dplyr::mutate(long = sf::st_coordinates(.)[,1], lat = sf::st_coordinates(.)[,2]) %>% st_drop_geometry() ################ PREPARE DATA FOR MODELLING ---------------------------------- DATA <- track_nofor_day_build %>% mutate(YDAY=yday(t_), hour=hour(t_), month=month(t_)) %>% filter(!is.na(step_length)) %>% filter(!is.na(turning_angle)) %>% filter(!is.na(speed)) %>% filter(!is.na(mean_speed)) %>% filter(!is.na(mean_angle)) %>% select(-tod_,-FOREST,-forest_size,-build_id,-type_of_food,-frequency) %>% mutate(point_id=seq_along(t_)) %>% mutate(year=as.numeric(year), bird_id=as.factor(bird_id)) %>% filter(!is.na(age_cy)) ##########~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###################################### ########## READ IN VALIDATION DATA FROM INDEPENDENT FEEDER SURVEYS ----------------- ##########~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~###################################### ## read in survey data from Eva Cereghetti ## Q1 is the question whether they feed or not #setwd("C:/Users/sop/OneDrive - Vogelwarte/REKI/Analysis/REKIfeeding") feed_surveys<-fread("data/survey.final.csv") %>% #filter(Q1=="Ja") %>% mutate(FEEDER_surveyed=ifelse(Q1=="Ja",1,0)) %>% mutate(ID=paste0("Eva",nr)) %>% select(ID,coord_x,coord_y,square,random,area,building.type,FEEDER_surveyed) %>% st_as_sf(coords = c("coord_x", "coord_y"), crs=21781) %>% st_transform(crs = 3035) ## read in survey from Fiona Pellet provided with addresses only ## Jerome Guelat provided R script to convert addresses to coordinates source("C:/Users/sop/OneDrive - Vogelwarte/General/ANALYSES/DataPrep/swisstopo_address_lookup.r") ## Feeding is the question whether they feed or not feed_surveys2<-read_csv("data/FeedersFionaPelle.csv", locale = locale(encoding = "UTF-8")) %>% #<-fread("data/FeedersFionaPelle.csv", encoding = 'UTF-8') %>% mutate(FEEDER_surveyed=ifelse(Feeding=="Yes",1,0)) ## generate coordinates from addresses feed_surveys2_locs<-swissgeocode(address=as.character(feed_surveys2$Address), nresults=3) feed_surv2_sf<-feed_surveys2_locs %>% rename(Address=address_origin) %>% left_join(feed_surveys2, by="Address",relationship ="many-to-many") %>% filter(!is.na(x)) %>% filter(!is.na(FEEDER_surveyed)) %>% mutate(ID=paste0("Fiona",ID)) %>% group_by(ID,lon,lat) %>% summarise(FEEDER_surveyed=max(FEEDER_surveyed)) %>% st_as_sf(coords = c("lon", "lat"), crs=4326) %>% st_transform(crs = 3035) %>% bind_rows(feed_surveys) %>% distinct() feed_surv2_sf$FEEDER_surveyed ## REMOVE DUPLICATE LOCATIONS mtx_distance <- st_distance(feed_surv2_sf[feed_surv2_sf$FEEDER_surveyed==1,], feed_surv2_sf[feed_surv2_sf$FEEDER_surveyed==1,]) mtx_distance<-as_tibble(mtx_distance) names(mtx_distance)<-feed_surv2_sf$ID[feed_surv2_sf$FEEDER_surveyed==1] eliminate<-mtx_distance %>% mutate(ID=feed_surv2_sf$ID[feed_surv2_sf$FEEDER_surveyed==1]) %>% gather(key="ID2", value="dist",-ID) %>% mutate(dist=as.numeric(dist)) %>% filter(!(ID==ID2)) %>% filter(dist<40) %>% arrange(desc(dist)) dim(feed_surv2_sf) for (l in 1:dim(eliminate)[1]) { #gridids<-length(unique(feed_surv2_sf$gridid[feed_surv2_sf$ID %in% c(eliminate$ID[l],eliminate$ID2[l])])) ## check whether the two points are in the same grid cell pair1<-eliminate$ID2[l] if((eliminate$ID[l] %in% feed_surv2_sf$ID)){ feed_surv2_sf<-feed_surv2_sf %>% filter(!(ID==pair1)) } } dim(feed_surv2_sf) export<-feed_surv2_sf %>% st_transform(4326) %>% dplyr::mutate(long = sf::st_coordinates(.)[,1], lat = sf::st_coordinates(.)[,2]) fwrite(export,"C:/Users/sop/OneDrive - Vogelwarte/General/MANUSCRIPTS/AnthropFeeding/DataArchive/REKI_validation_feeders.csv")