Skip to content

Instantly share code, notes, and snippets.

@jamesgrecian

jamesgrecian/kernel.r

Created July 17, 2018 13:39
Show Gist options
  • Save jamesgrecian/6a2f34ab85a8f9fbcca486997aebd0c8 to your computer and use it in GitHub Desktop.
Save jamesgrecian/6a2f34ab85a8f9fbcca486997aebd0c8 to your computer and use it in GitHub Desktop.
Download ZIP
Example script to generate a bivariate kernel from telemetry data using adehabitat
Raw
kernel.r
#####################################################################################################
### Example script to generate and plot a bivariate kernel utilisation distribution from GPS data ###
#####################################################################################################
# Load libraries
require(tidyverse)
require(viridis)
require(rworldmap)
require(sf)
require(raster)
require(adehabitatHR)
# Define projection - UTM 30N in this case
prj = "+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Define bounding box
CP <- as(raster::extent(-20, 15, 45, 65), "SpatialPolygons")
proj4string(CP) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
# Load in the world shapefile from rworldmap, clip to bounding box and project
proj4string(countriesLow) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
world_shp <- rgeos::gIntersection(countriesLow, CP, byid=T, drop_lower_td = T)
world_utm = sp::spTransform(world_shp, CRS(prj))
#Load in sample GPS data
dat <- structure(list(BTO = structure(c(12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L, 12L,
12L, 12L, 12L, 12L, 12L), .Label = c("1446262", "1446266", "1446268",
"1446269", "1446273", "1459917", "1459918", "1459919", "1459920",
"1459927", "1459929", "1459930", "1459933", "1459936", "1459937",
"1459938", "1459939", "1484901", "1484903", "1484904", "1484905",
"1484906", "1484911", "1484912", "1484918", "1484920", "1484921",
"1484922", "1484923", "1484924", "1484932"), class = "factor"),
DateTime = structure(c(1435648800, 1435649400, 1435650000,
1435650600, 1435651200, 1435651800, 1435652400, 1435653000,
1435653600, 1435654200, 1435654800, 1435655400, 1435656000,
1435656600, 1435657200, 1435657800, 1435658400, 1435659000,
1435659600, 1435660200, 1435660800, 1435661400, 1435662000,
1435662600, 1435663200, 1435663800, 1435664400, 1435665000,
1435665600, 1435666200, 1435666800, 1435667400, 1435668000,
1435668600, 1435669200, 1435669800, 1435670400, 1435671000,
1435671600, 1435672200, 1435672800, 1435673400, 1435674000,
1435674600, 1435675200, 1435675800, 1435676400, 1435677000,
1435677600, 1435678200, 1435678800, 1435679400, 1435680000,
1435680600, 1435681200, 1435681800, 1435682400, 1435683000,
1435683600, 1435684200, 1435684800, 1435685400, 1435686000,
1435686600, 1435687200, 1435687800, 1435688400, 1435689000,
1435689600, 1435690200, 1435690800, 1435691400, 1435692000,
1435692600, 1435693200, 1435693800, 1435694400, 1435695000,
1435695600, 1435696200, 1435696800, 1435697400, 1435698000,
1435698600, 1435699200, 1435699800, 1435700400, 1435701000,
1435701600, 1435702200, 1435702800, 1435703400, 1435704000,
1435704600, 1435705200, 1435705800, 1435706400, 1435707000,
1435707600, 1435708200, 1435708800, 1435709400, 1435710000,
1435710600, 1435711200, 1435711800, 1435712400, 1435713000,
1435713600, 1435714200, 1435714800, 1435715400, 1435716000,
1435716600, 1435717200, 1435717800, 1435718400, 1435719000,
1435719600, 1435720200, 1435720800, 1435721400, 1435722000,
1435722600, 1435723200, 1435723800, 1435724400, 1435725000,
1435725600, 1435726200, 1435726800, 1435727400, 1435728000,
1435728600, 1435729200, 1435729800, 1435730400, 1435731000,
1435731600, 1435732200, 1435732800, 1435733400, 1435734000,
1435734600, 1435735200, 1435735800, 1435736400, 1435737000,
1435737600, 1435738200, 1435738800, 1435739400, 1435740000,
1435740600, 1435741200, 1435741800, 1435742400, 1435743000,
1435743600, 1435744200, 1435744800, 1435745400, 1435746000,
1435746600, 1435747200, 1435747800, 1435748400, 1435749000,
1435749600, 1435750200, 1435750800, 1435751400, 1435752000,
1435752600, 1435753200, 1435753800, 1435754400, 1435755000,
1435755600, 1435756200), class = c("POSIXct", "POSIXt"), tzone = "GMT"),
lon = structure(c(-2.64632855434523, -2.59229330597376, -2.52932700577752,
-2.53029230689389, -2.47260001348875, -2.52230475745883,
-2.53397883030807, -2.4983039239302, -2.42944343248492, -2.32210523073143,
-2.21162796934958, -2.07619435183371, -1.94369668682915,
-1.80499868078535, -1.70662909587723, -1.69955568532313,
-1.57655531412398, -1.55864745818821, -1.54369867081544,
-1.53041896617127, -1.39176187113321, -1.2559091729464, -1.12045745124851,
-0.974796932798973, -0.828562489837168, -0.683781394796058,
-0.537671388483779, -0.395028637030832, -0.263843862246226,
-0.120806930634551, 0.0127831811227706, 0.143313715433408,
0.277003846682696, 0.411918480020646, 0.529367058447774,
0.672513481239938, 0.801679209471577, 0.922985847546624,
0.983010976443613, 1.08384311586926, 1.17291441418054, 1.21468100538223,
1.19911796994562, 1.24196338444466, 1.29846807105857, 1.31664045473779,
1.36732952017921, 1.36952270541004, 1.3642590071282, 1.36852771897301,
1.36797883524332, 1.36810502126107, 1.36432290264499, 1.33884534453655,
1.34297126980888, 1.32512457479391, 1.33565796603004, 1.34072907938896,
1.35275724184489, 1.34708271594895, 1.34701629716973, 1.34738399926366,
1.34783878283474, 1.35015685720195, 1.35377868355188, 1.35629913956653,
1.32093290065812, 1.30058451171807, 1.33099257716498, 1.33011248468776,
1.31027388678288, 1.30626461086156, 1.31244687331144, 1.31364654496076,
1.32826214257485, 1.34486535310086, 1.34310721940544, 1.36717566356309,
1.36091670727838, 1.3619685540277, 1.36184035681361, 1.2927789328389,
1.29848879947131, 1.2995681719765, 1.30076346972821, 1.30020647589823,
1.30060555850547, 1.30053331606335, 1.30138051322503, 1.30110653808958,
1.30054093933714, 1.30062144096979, 1.30069467912499, 1.3002894553097,
1.30048828124674, 1.29949993697569, 1.29854655978769, 1.29715665591914,
1.29899471623905, 1.29828456631371, 1.29471019298013, 1.29418674813264,
1.29258207243124, 1.29128872662572, 1.28978116495392, 1.28884108166261,
1.28804248590252, 1.28656937806221, 1.2854450127562, 1.28430821108267,
1.28366606601372, 1.28336330747078, 1.28325180767614, 1.28323872102095,
1.28331621582762, 1.28345585796003, 1.28348645134456, 1.2828814789377,
1.28251307830534, 1.28188708572819, 1.28121041691473, 1.28411418635269,
1.28388163679546, 1.27521254248526, 1.2728954083403, 1.27199422236911,
1.27270546438984, 1.27376861744277, 1.27398975338933, 1.27478315223322,
1.29226871909555, 1.29328387953101, 1.29910087261316, 1.34362869776121,
1.34005078334307, 1.34227248831918, 1.30763254003217, 1.29376285596822,
1.28197051759926, 1.25735350217047, 1.24476386640919, 1.25684114947309,
1.19939723968055, 1.19200237845431, 1.19200099598189, 1.18079527744403,
1.18780135925326, 1.20139207959947, 1.20513842186178, 1.21637829991362,
1.22018102807916, 1.14600067874507, 1.08908586042547, 1.00255949107928,
0.872015304750368, 0.741641339226441, 0.606138530064448,
0.471172719138869, 0.449838406957629, 0.384442162629911,
0.257016339973394, 0.123438313559227, -0.0368993395593446,
-0.212657344524883, -0.37361682210631, -0.529314013156887,
-0.678831365707043, -0.841852459773389, -1.00245590573399,
-1.10057402405274, -1.23233071078879, -1.38692176534051,
-1.54813853163459, -1.6769929275925, -1.79469545467347, -1.93378181668947,
-2.1062178546785, -2.27976118668737, -2.44461424592071, -2.61058763573889
), .Names = c("1", "2", "3", "4", "5", "6", "7", "8", "9",
"10", "11", "12", "13", "14", "15", "16", "17", "18", "19",
"20", "21", "22", "23", "24", "25", "26", "27", "28", "29",
"30", "31", "32", "33", "34", "35", "36", "37", "38", "39",
"40", "41", "42", "43", "44", "45", "46", "47", "48", "49",
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59",
"60", "61", "62", "63", "64", "65", "66", "67", "68", "69",
"70", "71", "72", "73", "74", "75", "76", "77", "78", "79",
"80", "81", "82", "83", "84", "85", "86", "87", "88", "89",
"90", "91", "92", "93", "94", "95", "96", "97", "98", "99",
"100", "101", "102", "103", "104", "105", "106", "107", "108",
"109", "110", "111", "112", "113", "114", "115", "116", "117",
"118", "119", "120", "121", "122", "123", "124", "125", "126",
"127", "128", "129", "130", "131", "132", "133", "134", "135",
"136", "137", "138", "139", "140", "141", "142", "143", "144",
"145", "146", "147", "148", "149", "150", "151", "152", "153",
"154", "155", "156", "157", "158", "159", "160", "161", "162",
"163", "164", "165", "166", "167", "168", "169", "170", "171",
"172", "173", "174", "175", "176", "177", "178", "179", "180"
)), lat = structure(c(56.0825894041409, 56.0864570290153,
56.1106450868391, 56.1119246164043, 56.127073939209, 56.1197377968335,
56.1165118881497, 56.1169292200395, 56.1529294297373, 56.200223056039,
56.2400094227424, 56.2589204112125, 56.3026783887496, 56.3420310398718,
56.3266578761411, 56.325793120515, 56.3333006680866, 56.3341629178721,
56.3279408224272, 56.3289072101196, 56.3255853500328, 56.3195506208492,
56.3176680619058, 56.3268477587496, 56.3361791758799, 56.3582265214522,
56.3748850527724, 56.4013053125873, 56.43966348614, 56.4828788130415,
56.524535722733, 56.5745171779971, 56.627716901841, 56.6730577956402,
56.7125229506047, 56.7458299870163, 56.7488319630706, 56.7955244419465,
56.8327732054423, 56.897685622386, 56.9418090963799, 56.9529519392587,
56.9578650242673, 57.0023379732139, 57.0460324805241, 57.0518820263515,
57.1190436107416, 57.1200733440476, 57.1221996826841, 57.1233009654003,
57.1262881591838, 57.1183150144299, 57.1187707228445, 57.1661369138986,
57.1662817497776, 57.165640042079, 57.17338840518, 57.1725652534226,
57.1732452021958, 57.1713016643192, 57.1756930066256, 57.1846471942915,
57.1873295934383, 57.1896145508359, 57.1933022109752, 57.1954001589519,
57.2077317403106, 57.1919550275635, 57.1679106197951, 57.1636329876302,
57.1692637182242, 57.1895925036955, 57.1935530206832, 57.1950586615873,
57.1964974429232, 57.2078758660872, 57.2065967012277, 57.2161305730806,
57.2073544663854, 57.2039271373198, 57.2048759667625, 57.2264344342153,
57.228334360384, 57.2297273496261, 57.2308195645804, 57.2301101783512,
57.2305461763228, 57.2311782102726, 57.2317877054901, 57.23211188055,
57.2322392012487, 57.232218453346, 57.2318651888746, 57.231789148031,
57.2311617672851, 57.2305930804566, 57.2298972766741, 57.229297368126,
57.2282360748448, 57.2278197288013, 57.2275554669875, 57.226790946341,
57.2261533770362, 57.2254691599883, 57.2248102462054, 57.2240689520879,
57.2235746343459, 57.2229684679149, 57.2223953348034, 57.2219575210611,
57.2220244278481, 57.2218728858636, 57.2218407203025, 57.2219850647048,
57.2220813469885, 57.2224273266866, 57.2229479096229, 57.2235351699995,
57.22437516172, 57.2250530700945, 57.2258945892989, 57.2312095578689,
57.2323057252657, 57.2337167455515, 57.2343443163611, 57.2357420314075,
57.2378582623771, 57.240133377364, 57.2416547421964, 57.2439480362934,
57.2350318169528, 57.2368137823532, 57.2386650767769, 57.2179019389071,
57.2101741944763, 57.207911744742, 57.1895185115749, 57.1837618621134,
57.1871062654934, 57.1999692222485, 57.2061069924648, 57.2074691889366,
57.1929939076002, 57.1919284388537, 57.1951570310224, 57.1966194461799,
57.2028426534667, 57.2129639133466, 57.2223081730905, 57.2353307057835,
57.2430539365511, 57.2395832012003, 57.2311277626122, 57.208128745782,
57.1641128205851, 57.126066131963, 57.0813109552427, 57.0484641420103,
57.035407208896, 57.0114464815411, 56.9594190013255, 56.9140724384926,
56.8796521369934, 56.8536979213514, 56.812718508252, 56.7733453901123,
56.7238706261454, 56.6862452561409, 56.6453512652781, 56.6121026794116,
56.5622237240516, 56.5086962209695, 56.4672049647153, 56.421526468584,
56.3702705472641, 56.3232681124185, 56.2669038622567, 56.2080162323506,
56.161507206751, 56.1017185026255), .Names = c("1", "2",
"3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13",
"14", "15", "16", "17", "18", "19", "20", "21", "22", "23",
"24", "25", "26", "27", "28", "29", "30", "31", "32", "33",
"34", "35", "36", "37", "38", "39", "40", "41", "42", "43",
"44", "45", "46", "47", "48", "49", "50", "51", "52", "53",
"54", "55", "56", "57", "58", "59", "60", "61", "62", "63",
"64", "65", "66", "67", "68", "69", "70", "71", "72", "73",
"74", "75", "76", "77", "78", "79", "80", "81", "82", "83",
"84", "85", "86", "87", "88", "89", "90", "91", "92", "93",
"94", "95", "96", "97", "98", "99", "100", "101", "102",
"103", "104", "105", "106", "107", "108", "109", "110", "111",
"112", "113", "114", "115", "116", "117", "118", "119", "120",
"121", "122", "123", "124", "125", "126", "127", "128", "129",
"130", "131", "132", "133", "134", "135", "136", "137", "138",
"139", "140", "141", "142", "143", "144", "145", "146", "147",
"148", "149", "150", "151", "152", "153", "154", "155", "156",
"157", "158", "159", "160", "161", "162", "163", "164", "165",
"166", "167", "168", "169", "170", "171", "172", "173", "174",
"175", "176", "177", "178", "179", "180")), states = structure(c(2L,
2L, 2L, 2L, 1L, 2L, 1L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 3L, 2L, 2L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L,
1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 2L,
2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 1L, 1L, 2L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L,
1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 2L, 3L, 3L, 3L, 3L, 3L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("1",
"2", "3"), class = "factor")), .Names = c("BTO", "DateTime",
"lon", "lat", "states"), row.names = c(NA, 180L), class = c("tbl_df",
"tbl", "data.frame"))
# The new sf package allows better manipulaton of spatial dataframes incl. improved plotting and on the fly transformations...
dat_sf <- st_as_sf(dat, coords = c("lon", "lat")) %>% st_set_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#Project to UTM
dat_sf <- st_transform(st_as_sf(dat_sf), prj)
#Use adehabitat function to generate kernel specifying the smoothing parameter h in the same units as the locations (metres in UTM)
UD = kernelUD(as_Spatial(dat_sf$geometry), h = 10000)
image(ud)
#Extract the vertices for nice shapefile plotting
ver25 = adehabitatHR::getverticeshr(UD, percent = 25, ida = "25%")
ver50 = adehabitatHR::getverticeshr(UD, percent = 50, ida = "50%")
ver75 = adehabitatHR::getverticeshr(UD, percent = 75, ida = "75%")
ver95 = adehabitatHR::getverticeshr(UD, percent = 95, ida = "95%")
# Then use spRbind to combine the polygons, can only combine two at a time
# Layer from the bottom up
a = maptools::spRbind(ver95, ver75)
b = maptools::spRbind(ver50, ver25)
output = maptools::spRbind(a, b)
#Generate a nice plot...
p1 <- ggplot() +
theme_bw() +
geom_sf(aes(fill = id, colour = id), data = st_as_sf(output)) +
geom_sf(data = st_transform(st_as_sf(world_shp), prj), fill = "dark grey") +
coord_sf(xlim = c(260000, 1250000), ylim = c(5750000, 6750000), crs = prj, expand = T) +
scale_fill_viridis_d() + scale_colour_viridis_d()
print(p1)
#ends
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment