Distribution of white browed sparrow weavers

Author

Martine van den Berg

Is the distribution of white-browed sparrow weavers (Plocepasser mahali) related to temperature or vegetation type?

To determine how the distribution of white-browed sparrow weavers is related to temperature and vegetation type, the following data will be used:

Observation data on specie distribution from Inaturalist
MODIS surface temperature
SANBI 2018 VegMap
South African Political border from ArcGIS Hub

The following R packages will be used for this analysis: sf(), terra(), rosm(), ggspatial(), tidyverse(), wesanderson()

library(sf)

Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(terra)

terra 1.7.3

library(rosm)
library(ggspatial)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.1     ✔ tibble    3.1.8
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks terra::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(wesanderson)

Plotting a map of white-browed sparrow weaver distribution and vegetation type

First, data on specie distribution from Inaturalist needs to be read in. This then needs to be be turned into a spatial object with the class sf, and the coordinate reference system (CRS) needs to be defined. Note that for this analysis I will use a geographic CRS since the entire South-Africa will be represented on the maps. This is a large area and projections such as Transverse Mercator would not be appropriate.

After reading in the data, the points are plotted on a basemap to get an initial look at the species distribution

#Read in SANBI VegMap and plot to look at data
obs <- read.csv("Data/inat_wbsw.csv")
#make into spatial object and define crs
wbsw <- st_as_sf(obs, coords = c("longitude", "latitude"), crs = 4326)

#Confirm that the data has the correct crs and class
st_crs(wbsw)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

class(wbsw)

[1] "sf"         "data.frame"

#What are the attribute columns?
names(wbsw)

 [1] "id"                               "observed_on_string"              
 [3] "observed_on"                      "time_observed_at"                
 [5] "time_zone"                        "user_id"                         
 [7] "user_login"                       "user_name"                       
 [9] "created_at"                       "updated_at"                      
[11] "quality_grade"                    "license"                         
[13] "url"                              "image_url"                       
[15] "sound_url"                        "tag_list"                        
[17] "description"                      "num_identification_agreements"   
[19] "num_identification_disagreements" "captive_cultivated"              
[21] "oauth_application_id"             "place_guess"                     
[23] "positional_accuracy"              "private_place_guess"             
[25] "private_latitude"                 "private_longitude"               
[27] "public_positional_accuracy"       "geoprivacy"                      
[29] "taxon_geoprivacy"                 "coordinates_obscured"            
[31] "positioning_method"               "positioning_device"              
[33] "species_guess"                    "scientific_name"                 
[35] "common_name"                      "iconic_taxon_name"               
[37] "taxon_id"                         "geometry"

#Plot points on a basemap
ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data=wbsw)

I only want the observations from South Africa, and will use the function st_intersection(), to achieve this., I use the South African Political border, which was downloaded from ArcGIS Hub.

#Read in data on SA political border
sa <- st_read("Data/sa_boundary/sa_poly.shp")

Reading layer `sa_poly' from data source 
  `C:\Users\Martine\Documents\R_work\GIS_deliverable\Data\sa_boundary\sa_poly.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 16.45306 ymin: -34.83427 xmax: 32.89441 ymax: -22.12656
Geodetic CRS:  WGS 84

# Total observations 
dim(wbsw)

[1] 989  38

#Create object of white browed sparrow weavers in SA only
wbsw2 <- st_intersection(wbsw, sa)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

dim(wbsw2)

[1] 462  39

#Observations after intersection
dim(wbsw)

[1] 989  38

#Plot
ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data=wbsw2)

Next, the SANBI VegMap data needs to be read into R. Then dissolve adjacent polygons with shared boundaries, so the VegMap can be plotted without the unwanted boundaries. I am interested in seeing the different biomes in which white-browed sparrow weavers occur, so the VegMap data will be grouped by biome.

#Read in SANBI VegMap 
vegmap<-st_read("Data/SANBI_VEGMAP2018/NVM2018_AEA_V22_7_16082019_final.shp")

Reading layer `NVM2018_AEA_V22_7_16082019_final' from data source 
  `C:\Users\Martine\Documents\R_work\GIS_deliverable\Data\SANBI_VEGMAP2018\NVM2018_AEA_V22_7_16082019_final.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 48732 features and 28 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -833551.4 ymin: -3766966 xmax: 781623.8 ymax: -2356961
Projected CRS: AEA_RSA_WGS84

#Dissolve
veg<-vegmap %>% group_by(BIOME_18) %>% 
  summarize()

#Confirm CRS
st_crs(veg)

Coordinate Reference System:
  User input: AEA_RSA_WGS84 
  wkt:
PROJCRS["AEA_RSA_WGS84",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Albers Equal Area",
            ID["EPSG",9822]],
        PARAMETER["Latitude of false origin",0,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",25,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",-24,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",-33,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

#Plot
veg %>% 
  ggplot() + geom_sf(aes(fill = BIOME_18))

Finally, I can create a map showing both the distribution of white browed sparrow weavers, and the different biomes in South Africa

#Create map of wbsw and vegmap
veg %>% 
  ggplot() + 
  geom_sf(aes(fill = BIOME_18))+
  geom_sf(data=wbsw2)

Creating a map of white-browed sparrow weaver distribution and surface temperature

First, I need to read in the MODIS surface temperature data, which was extracted with Java Script from Google Earth Engine. Then look at the CRS and class

#Read in temp data and look at data
temp <- rast("Data/Modis-Longterm-mean-LST-day.tif")

class(temp)

[1] "SpatRaster"
attr(,"package")
[1] "terra"

temp

class       : SpatRaster 
dimensions  : 1415, 1831, 1  (nrow, ncol, nlyr)
resolution  : 0.008983153, 0.008983153  (x, y)
extent      : 16.44815, 32.89631, -34.83667, -22.12551  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : Modis-Longterm-mean-LST-day.tif 
name        : LST_Day_1km 
min value   :    17.89276 
max value   :    38.09803

#What is the CRS?
crs(temp)

[1] "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"

Finally, I can create the map of white-browed sparrow weavers along with the surface temperature data for South-Africa

#Create an object of the colour palette
pal<-wes_palette("Zissou1", 100, type = "continuous")

#Plot
tempplot<-ggplot() +
    geom_raster(data = as.data.frame(temp, xy = TRUE),
                aes(x = x, y = y,
                    fill=LST_Day_1km))+
                    scale_fill_gradientn(colours = pal)+
  geom_sf(data=wbsw2, colour = "black", size = 1)+
  labs(x="Latitude",y="Longitude")

#Change label of key
tempplot+labs(fill = "Surface Temperature ˚C")

References for data used

Esri (2023). “Feature Layer” [basemap]. Scale 1:5,000,000. “World Countries (Generalized)”. Accessed 2023-03-10 from https://hub.arcgis.com/datasets/esri::world-countries-generalized

iNaturalist community. Observations of Plocepasser mahali from southern Africa observed between 2006-2022. Accessed 2023-03-13 from https://www.inaturalist.org

South African National Biodiversity Institute (2006). The Vegetation Map of South Africa, Lesotho and Swaziland, Mucina, L., Rutherford, M.C. and Powrie, L.W. (Editors), Online, Version 2018 http://bgis.sanbi.org/SpatialDataset/Detail/18,

Wan, Z., Hook, S., Hulley, G., (2021) MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid VO61 [Data set] NASA EOSDIS Land Processes DAAC. Accessed 2023-03-10 from https://doi.org/10.5067/MODIS/MOD11A1.061