knitr::opts_chunk$set()
# load libraries
library(raster)
library(rgdal)
library(spatstat)
library(maptools)
library(tidyverse)
library(mapview)
#library(qgisprocess)
memory.limit(size=56000)
## [1] 56000
#options(scipen=999)

Bulinus

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R\ShireRiver", layer: "bulinus"
## with 17 features
## It has 13 fields
## Integer64 fields read as strings:  field_1 Bulinus_pr Biomphalar
{plot(temperature,
      ylab = 'UTM Zone 36 S Northing [m]', 
      xlab = 'UTM Zone 36 S Easting [m]', 
      col = pal(6))
plot(bulinus, add = TRUE, pch = "+")}

elevation <- raster("Malawi_SRTM30meters.tif") 
  
elevation <- crop(elevation, extent(southern_mw)) |> 
  mask(southern_mw)
# project raster
raster::crs(elevation) <- proj4string(temperature)
#crs(elevation) <- "+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs"
# compareRaster(elevation, temperature) FALSE
extent(elevation) <- extent(temperature)
elevation_img<- as.im(elevation)

{plot(elevation,
      ylab = 'UTM Zone 36 S Northing [m]', 
      xlab = 'UTM Zone 36 S Easting [m]', 
      col = rev(terrain.colors(255)))
plot(bulinus, add = TRUE, pch = "+")}

mapView(elevation,
        trim = TRUE)+
mapview(bulinus, zcol= "Bulinus_co",
        layer.name = "Bulinus")
# function to convert spatial points to planar point pattern
convert.ppp <- function(points){
  x <- coordinates(points)[,1]
  y <- coordinates(points)[,2]
  l1 <- bbox(points)[1,1]
  l2 <- bbox(points)[1,2]
  l3 <- bbox(points)[2,1]
  l4 <- bbox(points)[2,2]
  win <- owin(xrange = c(l1,l2), yrange = c(l3,l4))
  points.ppp <- spatstat.geom::ppp(x, y,
                                   window = win,
                                   marks = points@data)
  return(points.ppp)
}

bulinus.ppp <- convert.ppp(bulinus)
unitname(bulinus.ppp)<- c("degree", "degrees")


# observation window from extent
window <- owin(xrange = c(temperature@extent[1:2]),
               yrange = c(temperature@extent[3:4]))
Window(bulinus.ppp) <- window
bulinus.ppp <- unmark(bulinus.ppp)

# identification of covariates 
isCovariate <- function(var, ppp){
  b <- quantile(var, probs = (0:4)/4)
  cut <- cut(var, breaks = b, label = 1:4) 
  tess <- tess(image = cut)  
  plot(tess)
  qq <- quadratcount(ppp, tess = tess)
  plot(qq, add = TRUE, pch = "+")
  quadrat.test(qq)
}

isCovariate(temperature, bulinus.ppp)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 1.5897, df = 3, p-value = 0.6765
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)
isCovariate(elevation, bulinus.ppp)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 29.122, df = 3, p-value = 4.222e-06
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)

It appears that: one) Bulinus. is abundant in moderately-elavated areas (at 100m and 500m, in rising order) and areas with mild-hot temperatures (particularly around 28-33 degrees celsius), two) temperature offers little explanation for the distribution of Bulinus i.e., seems it is a weak covariate, based on the pilot survey data, and three) elevation is a strong influencer of Bulinus abundance.

# covariate-based prediction
covariatePrediction <- function(var, ppp, label){
  rho <- rhohat(ppp, var)
  plot(rho,
       xlab = label)
  pred <- predict(rho)
  plot(pred)
}

covariatePrediction(temperature_img, bulinus.ppp, 
                    label = "Temperature (Celsius)")

covariatePrediction(elevation_img, bulinus.ppp,
                    label = "Elevation (m)")

Biomphalaria

# load Biomphalaria locations
biomphalaria <- readOGR(dsn = getwd(), layer = "biomphalaria")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\R\ShireRiver", layer: "biomphalaria"
## with 6 features
## It has 13 fields
## Integer64 fields read as strings:  field_1 Bulinus_pr Biomphalar
{plot(temperature,
      ylab = 'UTM Zone 36 S Northing [m]', 
      xlab = 'UTM Zone 36 S Easting [m]', 
      col = pal(6))
plot(biomphalaria, add = TRUE, pch = "+")}

{plot(elevation,
      ylab = 'UTM Zone 36 S Northing [m]', 
      xlab = 'UTM Zone 36 S Easting [m]', 
      col = rev(terrain.colors(255)))
plot(biomphalaria, add = TRUE, pch = "+")}

biomphalaria.ppp <- convert.ppp(biomphalaria)
unitname(biomphalaria.ppp)= c("degree", "degrees")


# observation window from extent
window <- owin(xrange = c(temperature@extent[1:2]),
               yrange = c(temperature@extent[3:4]))
Window(biomphalaria.ppp) <- window
biomphalaria.ppp <- unmark(biomphalaria.ppp)


isCovariate(temperature, biomphalaria.ppp)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 3.3338, df = 3, p-value = 0.6859
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)
isCovariate(elevation, biomphalaria.ppp)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 16.902, df = 3, p-value = 0.001481
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)
covariatePrediction(temperature_img, biomphalaria.ppp, 
                    label = "Temperature (Celsius)")

covariatePrediction(elevation_img, biomphalaria.ppp,
                    label = "Elevation (m)")