# ==============================================================================
# * Libraries
# ==============================================================================
library(leaflet)
library(raster)
library(htmlwidgets)
library(dplyr)
library(sf)
library(sp)
library(rgdal)
library(move)
library(rgeos)
library(data.table)
library(adehabitatHR)
library(DT)
library(tidyr)
source("R/00_fxns.R")
eval_again = FALSE # setting to FALSE for lines that have been executed already
# and which I don't want to run each time I knit# ==============================================================================
# * Read in GPS data with foraging info
# ==============================================================================
eobs = readRDS("data/gps_forage.RDS") %>%
rbindlist() %>%
filter(forage_pts) %>%
as.data.frame()
colnames(eobs) [ colnames(eobs) %in% "location.lat"] = "lat"
colnames(eobs) [ colnames(eobs) %in% "location.long"] = "lon"
lotek = readRDS("data/gps_lotek.RDS")
# lotek$forage_pts = NA
colnames(lotek) [ colnames(lotek) %in% "Latitude"] = "lat"
colnames(lotek) [ colnames(lotek) %in% "Longitude"] = "lon"
cols = names(lotek)[ names(lotek) %in% names(eobs)]
lotek = lotek[ , cols]
eobs = eobs[ , cols]
gps_pts = rbind(lotek[ , cols], eobs[ , cols])
# note, lotek has no forage pts, hence they're all NA
p = st_as_sf(gps_pts, coords = c("lon", "lat"), crs=4326)
p = as(p, "Spatial")
if(eval_again) {
# Export shapefile for GRASS to work with
writeOGR(p, dsn="data/gis/gps_pts", layer="gps_pts",
driver="ESRI Shapefile", overwrite_layer = TRUE)
}
# ==============================================================================
# * WDPA
# ==============================================================================
wdpa = readOGR("data/gis/WDPA_Feb2019_TZA-shapefile/",
layer = "WDPA_Feb2019_TZA-shapefile-polygons")
wdpa = wdpa[ wdpa$MARINE == 0, ]
cols_keep = c("NAME","DESIG_ENG","DESIG_TYPE","IUCN_CAT",
"STATUS", "STATUS_YR", "GOV_TYPE", "OWN_TYPE", "MANG_AUTH",
"VERIF","METADATAID")
wdpa = wdpa[ cols_keep ]
str(wdpa@data) # all are factors
wdpa@data[] = lapply(wdpa@data, as.character)
wdpa$STATUS_YR = as.numeric(wdpa$STATUS_YR)
map_wdpa = function(poly) {
map = leaflet() %>% addTiles() %>%
addPolygons(data = poly,
weight = 0.2,
popup = ~paste0("name: ", poly$NAME, "<br>",
"desig: ", wdpa$DESIG_ENG, "<br>",
"MANAG_AUTH: ", poly$MANG_AUTH))
map
}
map_wdpa(wdpa)
# * Subset wdpa
wdpa_export = wdpa[ wdpa$DESIG_ENG %in% c("Game reserve",
"Nature Reserve",
"Forest Reserve",
"National Park"), ]
wdpa_export = wdpa_export[ wdpa_export$MANG_AUTH != "Not Reported", ]
map_wdpa(wdpa_export)
# Export WDPA shapefile for GRASS to work with
if(eval_again) {
writeOGR(wdpa_export, dsn="data/gis/WDPA/", layer="wdpa_prune",
driver="ESRI Shapefile")
}
# proj string: +proj=utm +zone=37 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs
# from: http://spatialreference.org/ref/epsg/wgs-84-utm-zone-37s/
dir="/Users/nistara/projects/papers/fruit-bat/data/gis/guf_tza"
cd $dir
file_wgs="12mCrp_tza.tif"
file_utm="tza_utm.tif"
file_100m="tza_utm_100m.tif"
gdalinfo $file_wgs
# Change the CRS from WGS84 to UTM 37S
gdalwarp -t_srs '+proj=utm +zone=37 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs' $file_wgs $file_utm
gdalinfo $file_utm
# Change resolution to 100 by 100 m
gdalwarp -tr 100 -100 -r max $file_utm $file_100m
gdalinfo $file_100m
# delete big $file_utm
rm $file_utm
# make polygons
# ------------------------------------------------------------------------------
# distances can be calculated from rasters as well, but I'm creating polygons
# and subsetting built-up areas for visualization, etc as well.
gdal_polygonize.py $file_100m tza_guf_poly.shp
# ==============================================================================
# * Getting distances of GUF and WDPA from GPS points
# ==============================================================================
# New grass location based on GUF layer (to match extent and projection)
# ==============================================================================
grass70 -c ~/projects/papers/fruit-bat/data/gis/guf_tza/tza_guf_poly.shp ~/grassdata/paper_fruit-bat/
# *** Open grass location created above
# ------------------------------------------------------------------------------
grass70 ~/grassdata/paper_fruit-bat/PERMANENT
# *** g.gisenv
# LOCATION_NAME=bat-track_tza
# GISDBASE=/Users/nistara/grassdata
# MAPSET=PERMANENT
# GUI=text
# PID=73234
# *** g.proj -p
# -PROJ_INFO-------------------------------------------------
# name : WGS_1984_UTM_Zone_37S
# datum : wgs84
# ellps : wgs84
# proj : utm
# zone : 37
# south : defined
# no_defs : defined
# -PROJ_UNITS------------------------------------------------
# unit : Meter
# units : Meters
# meters : 1
# *** Pre-GRASS setup
# ------------------------------------------------------------------------------
# Refer to src/project_guf.sh
# 1. It changes CRS from WGS84 to UTM 37S.
# out = tza_utm.tif
# 2. Change resoution to 100 by 100 m (from 12 by 12 m)
# out = tza_utm_100m.tif
# *** Import guf_tza_poly_utm37S
# ------------------------------------------------------------------------------
v.in.ogr input=/Users/nistara/projects/papers/fruit-bat/data/gis/guf_tza/tza_guf_poly.shp layer=tza_guf_poly output=guf_tza_poly_utm37S
# *** Import gps_pts and WDPA (pruned)
# ------------------------------------------------------------------------------
# gps_pts is combined e-obs and lotek GPS fixes
v.import input=/Users/nistara/projects/papers/fruit-bat/data/gis/gps_pts/gps_pts.shp --overwrite
v.import input=/Users/nistara/projects/papers/fruit-bat/data/gis/WDPA/wdpa_prune.shp --overwrite
g.list vect -p
# *** Extract urban areas from guf polygon
# ------------------------------------------------------------------------------
v.db.select guf_tza_poly_utm37S | head
# NOTE: DN is the column name with the values (urban = 255)
v.extract input=guf_tza_poly_utm37S output=guf_urban -d where="(DN = 255)" --overwrite
# *** Get distances
# ------------------------------------------------------------------------------
v.db.select gps_pts | head
v.db.select guf_urban | head
v.db.select wdpa_prune | head
v.db.addcolumn gps_pts columns="dist_guf double precision"
v.distance from=gps_pts to=guf_urban upload=dist column=dist_guf
v.db.addcolumn gps_pts columns="dist_wdpa double precision","wdpa_name text"
v.distance from=gps_pts to=wdpa_prune upload=dist,to_attr \
column=dist_wdpa,wdpa_name to_colum=NAME
# *** Export csv file with distances
# ------------------------------------------------------------------------------
csv_dir=~/projects/papers/fruit-bat/data/
rm $csv_dir/dist_gps.csv
db.out.ogr input=gps_pts output=$csv_dir/dist_gps.csv \
format=CSV
# ==============================================================================
# * Read in data with distance info (done in GRASS)
# ==============================================================================
gps_dist = read.csv("data/dist_gps.csv", stringsAsFactors = FALSE)
moro_roost = which(gps_dist$roost)
df = gps_dist[ -moro_roost, ]
df$site = ifelse(grepl("M", df$tagID), "Morogoro", "Kilombero")
# kilombero
kilo = df[ df$site %in% "Kilombero", ]
kilo_eobs = kilo[ !is.na(kilo$forage_pts), ]
kilo_lotek = kilo[ is.na(kilo$forage_pts), ]
# morogoro
moro = df[ df$site %in% "Morogoro", ]
dist_l = list(kilo_eobs = kilo_eobs, kilo_lotek = kilo_lotek, moro = moro)
# ** Information about forage/feed-roost fixes in built-ip areas
# ------------------------------------------------------------------------------
guf_dist = lapply(dist_l, function(df) {
data.frame(dist_0 = sum(df$dist_guf == 0),
dist_100 = sum(df$dist_guf <= 100),
all_forage = nrow(df),
percent_100 = sum(df$dist_guf <= 100)/nrow(df) * 100)
})
mapply(cbind, guf_dist, "site" = names(dist_l), SIMPLIFY=F) %>% rbindlist## dist_0 dist_100 all_forage percent_100 site
## 1: 126 232 505 45.94059 kilo_eobs
## 2: 27 30 34 88.23529 kilo_lotek
## 3: 57 77 105 73.33333 moroguf_ind_bat = lapply(dist_l, function(df) {
df_summary = df %>%
dplyr::select(tagID, dist_guf) %>%
group_by(tagID) %>%
summarize(dist_100 = sum(dist_guf <= 100),
all_forage = n(),
percent_100 = sum(dist_guf <= 100)/all_forage * 100) %>%
as.data.frame()
})
guf_ind_bat## $kilo_eobs
## tagID dist_100 all_forage percent_100
## 1 K5309 40 150 26.66667
## 2 K5310 10 84 11.90476
## 3 K5311 21 25 84.00000
## 4 K5312 35 45 77.77778
## 5 K5313 99 108 91.66667
## 6 K5315 0 12 0.00000
## 7 K5317 0 48 0.00000
## 8 K5319 27 33 81.81818
##
## $kilo_lotek
## tagID dist_100 all_forage percent_100
## 1 K166359 10 11 90.90909
## 2 K166361 8 8 100.00000
## 3 K166363 0 3 0.00000
## 4 K166366 12 12 100.00000
##
## $moro
## tagID dist_100 all_forage percent_100
## 1 M166358 19 25 76.00000
## 2 M166362 2 4 50.00000
## 3 M166364 4 4 100.00000
## 4 M166365 6 6 100.00000
## 5 M166367 7 9 77.77778
## 6 M166369 11 11 100.00000
## 7 M166370 15 15 100.00000
## 8 M166371 11 11 100.00000
## 9 M166372 2 20 10.00000# ** Information about forage/feed-roost fixes in protected areas
# ------------------------------------------------------------------------------
wdpa_dist = lapply(dist_l, function(df) {
wdpa_df = df[ df$dist_wdpa <= 100, ]
df_n = nrow(df)
wdpa_df %>%
group_by(wdpa_name) %>%
summarize(n= n()) %>%
as.data.frame() %>%
mutate(
all_forage = df_n,
percent_n = n/df_n * 100)
})
mapply(cbind, wdpa_dist[ c(1, 3)], "site" = names(dist_l)[c(1, 3)], SIMPLIFY=F) %>%
rbindlist()## wdpa_name n all_forage percent_n site
## 1: Mikumi National Park 12 505 2.376238 kilo_eobs
## 2: Udzungwa Mountains National Park 152 505 30.099010 kilo_eobs
## 3: Nyandira 11 105 10.476190 morowdpa_areas = lapply(dist_l, function(df) {
df_summary = df %>%
group_by(tagID) %>%
summarize(total = n())
wdpa_df = df[ df$dist_wdpa <= 100, ]
df_n = nrow(df)
wdpa_df %>%
group_by(tagID, wdpa_name) %>%
summarize(n= n()) %>%
as.data.frame() %>%
left_join(df_summary, "tagID") %>%
mutate(percent_n = n/total * 100)
})
wdpa_areas## $kilo_eobs
## tagID wdpa_name n total percent_n
## 1 K5309 Udzungwa Mountains National Park 81 150 54.00000
## 2 K5310 Udzungwa Mountains National Park 58 84 69.04762
## 3 K5312 Udzungwa Mountains National Park 8 45 17.77778
## 4 K5313 Udzungwa Mountains National Park 5 108 4.62963
## 5 K5317 Mikumi National Park 12 48 25.00000
##
## $kilo_lotek
## [1] tagID wdpa_name n total percent_n
## <0 rows> (or 0-length row.names)
##
## $moro
## tagID wdpa_name n total percent_n
## 1 M166358 Nyandira 9 25 36
## 2 M166372 Nyandira 2 20 10# NOTE: foraging is near urban, but it also means that there's s much fragmentation that urban areas are near protected areas