Mapping Marajo Houses for NDVI

Set Up

rm(list = ls())
setwd("C:/Users/katie/OneDrive/Marajo Project")

x <- c("sp", "rgdal", "rgeos", "maptools", "sp", "sf", "raster", "MODIS", "gdalUtils", 
    "ggplot2", "RColorBrewer", "ggmap", "mapview", "openxlsx", "leaflet", "tidyverse", 
    "mapedit")


lapply(x, library, character.only = TRUE)

register_google(key = "AIzaSyDnw8wZpgPWmXmD4r3qTv6vkL0-Kx-3THk")

Read in House x y Data

P4S <- CRS("+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

setwd("C:/Users/katie/OneDrive/Marajo Project")

HouseDataRaw <- read.csv("Marajo House Data.csv")
HouseDataRaw <- subset(HouseDataRaw, HouseDataRaw$longitude != 0 & HouseDataRaw$latitude != 0)
colnames(HouseDataRaw)[colnames(HouseDataRaw) %in% c("longitude", "latitude")] <- c("lon", "lat")
HouseDataRaw <- select(HouseDataRaw, group, lon, lat)



HD <- SpatialPointsDataFrame(coords = cbind(HouseDataRaw$lon, HouseDataRaw$lat), 
                             proj4string = P4S, data = HouseDataRaw)
HD.sp <- SpatialPoints(coords = cbind(HouseDataRaw$lon, HouseDataRaw$lat), 
                       proj4string = P4S)
HD.sf <- st_as_sf(HouseDataRaw, coords = c("lon", "lat"), crs =  4269) 

Basemap using getData

BRmap <- getData('GADM', country = 'BR', download = TRUE, path = "/Users/katie/OneDrive/Marajo Project/" , level=0)



BR <- readRDS(paste0('/Users/katie/OneDrive/Marajo Project/', "gadm36_BRA_0_sp.rds"))
BR.sp<- gSimplify(BR, tol = 0.001, topologyPreserve = TRUE)
BR.sp <- SpatialPolygonsDataFrame(BR.sp, data = BR@data)
BR.sf <- as(BR.sp, "sf")


ggplot() + 
  geom_sf(data = HD.sf, fill = "light grey")+
  geom_sf(data = BR.sf)+
  theme_minimal()

# getData only gave map of whole country. I tried to find shapefile for just Marajo couldn't so have tried some google maps based approaches

Basemap using get_map

BR_basemap <- get_map(location=c(lon = -48.59373 , lat =-0.8294666 ), zoom= 11, maptype = 'satellite', source = 'google')



ggmap(BR_basemap)+
  geom_point(data = HouseDataRaw, aes(x = lon, y = lat, alpha = 0, colour = group)) +
  scale_x_continuous(limits = c(-48.7, -48.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1, -0.7), expand = c(0, 0))

Basemap using interactive leaflet and mapview

mapview(HD.sf, burst = TRUE)