Reading in the data:
#POINTS & LINES
#Reading in data
track_lines <- st_read("C:/Users/carme/NR6950/Final Project Part 3/data/Elephant_Research_Collar_09105/lines.shp")
## Reading layer `lines' from data source
## `C:\Users\carme\NR6950\Final Project Part 3\data\Elephant_Research_Collar_09105\lines.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 12.357 ymin: 8.001 xmax: 14.366 ymax: 8.549
## Geodetic CRS: WGS 84
#checking CRS
st_crs(track_lines)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
#reprojecting to 33N, as communicated by the file's metadata
track_utm <- st_transform(track_lines, crs = 32633)
st_crs(track_utm)
## Coordinate Reference System:
## User input: EPSG:32633
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["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)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Navigation and medium accuracy spatial referencing."],
## AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]
#Reading in data
track_points <- st_read("C:/Users/carme/NR6950/Final Project Part 3/data/Elephant_Research_Collar_09105/points.shp")
## Reading layer `points' from data source
## `C:\Users\carme\NR6950\Final Project Part 3\data\Elephant_Research_Collar_09105\points.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 911 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 12.357 ymin: 8.001 xmax: 14.366 ymax: 8.549
## Geodetic CRS: WGS 84
#checking CRS
st_crs(track_points)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
#reprojecting to 33N, as communicated by the file's metadata
track_pts_utm <- st_transform(track_points, crs = 32633)
st_crs(track_pts_utm)
## Coordinate Reference System:
## User input: EPSG:32633
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["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)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Navigation and medium accuracy spatial referencing."],
## AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]
plot(track_pts_utm['timestamp'], xlab = "Latitude", ylab = "Longitude", main = "Elephant Tracking Point Data, by Timestamp")
#LANDSAT
#Band 4 is Landsat 7's near-infrared band.
nir_jan <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020109_20200917_02_T1/LE07_L1TP_184054_20020109_20200917_02_T1_B4.TIF")
nir_mar <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020314_20200916_02_T1/LE07_L1TP_184054_20020314_20200916_02_T1_B4.TIF")
nir_may <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020517_20200916_02_T1/LE07_L1TP_184054_20020517_20200916_02_T1_B4.TIF")
nir_jul <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020720_20200916_02_T1/LE07_L1TP_184054_20020720_20200916_02_T1_B4.TIF")
nir_sep <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020922_20200916_02_T1/LE07_L1TP_184054_20020922_20200916_02_T1_B4.TIF")
nir_oct <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20021008_20200916_02_T1/LE07_L1TP_184054_20021008_20200916_02_T1_B4.TIF")
nir_nov <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20021125_20200916_02_T1/LE07_L1TP_184054_20021125_20200916_02_T1_B4.TIF")
#Band 4 is Landsat 7's red band.
red_jan <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020109_20200917_02_T1/LE07_L1TP_184054_20020109_20200917_02_T1_B3.TIF")
red_mar <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020314_20200916_02_T1/LE07_L1TP_184054_20020314_20200916_02_T1_B3.TIF")
red_may <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020517_20200916_02_T1/LE07_L1TP_184054_20020517_20200916_02_T1_B3.TIF")
red_jul <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020720_20200916_02_T1/LE07_L1TP_184054_20020720_20200916_02_T1_B3.TIF")
red_sep <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20020922_20200916_02_T1/LE07_L1TP_184054_20020922_20200916_02_T1_B3.TIF")
red_oct <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20021008_20200916_02_T1/LE07_L1TP_184054_20021008_20200916_02_T1_B3.TIF")
red_nov <- rast("C:/Users/carme/NR6950/Final Project Part 3/data/NDVI_data/LE07_L1TP_184054_20021125_20200916_02_T1/LE07_L1TP_184054_20021125_20200916_02_T1_B3.TIF")
#BENOUE AOI
benoue_np <- vect("C:/Users/carme/NR6950/Final Project Part 3/Benoue/Benoue_NP/Benoue.shp")
## Warning: [vect] Z coordinates ignored
benoue_33n <- project(benoue_np, crs(nir_jan))
Preparing the data:
#POINTS:
#Grab 2002 points
points_2002 <- track_pts_utm[year(track_pts_utm$timestamp) == 2002, ]
#Converting timestamp field to date and time
points_2002$timestamp <- as.POSIXct(points_2002$timestamp)
points_2002 <- points_2002 %>%
mutate(month = factor(month(timestamp, label = TRUE, abbr = TRUE),
levels = c("Jan","Feb","Mar","Apr","May","Jun","Jul",
"Aug","Sep","Oct","Nov","Dec")))
#Note: track_lines is a LINESTRING and does not have timestamp data.
#Cropping points to Benoue NP shapefile I made
points_sf <- st_as_sf(points_2002,
coords = c("lon", "lat"),
crs = 32633)
benoue_sf <- st_as_sf(benoue_33n)
points_sf <- st_transform(points_sf, st_crs(benoue_sf))
st_crs(benoue_sf)
## Coordinate Reference System:
## User input: WGS 84 / UTM zone 33N
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Navigation and medium accuracy spatial referencing."],
## AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]
st_crs(points_sf)
## Coordinate Reference System:
## User input: WGS 84 / UTM zone 33N
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Navigation and medium accuracy spatial referencing."],
## AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]
points <- st_intersection(points_sf, benoue_sf)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
#LANDSAT:
#Listing NIR rasters
nir_list <- list(nir_jan, nir_mar, nir_may, nir_jul, nir_sep, nir_oct, nir_nov)
exts <- lapply(nir_list, ext)
#extents
xmin <- max(sapply(exts, function(e) e[1]))
xmax <- min(sapply(exts, function(e) e[2]))
ymin <- max(sapply(exts, function(e) e[3]))
ymax <- min(sapply(exts, function(e) e[4]))
if (xmin >= xmax || ymin >= ymax) stop("No overlapping area among rasters!")
common_ext <- ext(xmin, xmax, ymin, ymax)
nir_list_cropped <- lapply(nir_list, function(r) crop(r, common_ext))
## |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|=========================================
sapply(nir_list_cropped, ext)
## [[1]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[2]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[3]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[4]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[5]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[6]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[7]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
nir_stack <- rast(nir_list_cropped)
names(nir_stack) <- c("Jan","Mar","May","Jul","Sep","Oct","Nov")
#Cropping NIR layers to Benoue NP shapefile I made
#cropping to Benoue NP polygon
nir_list_cropped_masked <- lapply(nir_list, function(r) {
r2 <- crop(r, benoue_33n)
mask(r2, benoue_33n)
})
sapply(nir_list_cropped_masked, function(x) c(ext(x), ncol(x), nrow(x)))
## [,1]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,2]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,3]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,4]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,5]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,6]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,7]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
#Stacking
nir_stack <- rast(nir_list_cropped_masked)
names(nir_stack) <- c("Jan","Mar","May","Jul","Sep","Oct","Nov")
#Listing red rasters
red_list <- list(red_jan, red_mar, red_may, red_jul, red_sep, red_oct, red_nov)
exts <- lapply(red_list, ext)
#extents
xmin <- max(sapply(exts, function(e) e[1]))
xmax <- min(sapply(exts, function(e) e[2]))
ymin <- max(sapply(exts, function(e) e[3]))
ymax <- min(sapply(exts, function(e) e[4]))
if (xmin >= xmax || ymin >= ymax) stop("No overlapping area among rasters!")
common_ext <- ext(xmin, xmax, ymin, ymax)
red_list_cropped <- lapply(red_list, function(r) crop(r, common_ext))
## |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|=========================================
sapply(red_list_cropped, ext)
## [[1]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[2]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[3]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[4]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[5]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[6]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
##
## [[7]]
## SpatExtent : 327585, 558015, 855585, 1062915 (xmin, xmax, ymin, ymax)
red_stack <- rast(red_list_cropped)
names(red_stack) <- c("Jan","Mar","May","Jul","Sep","Oct","Nov")
#Cropping red layers to Benoue NP shapefile I made
benoue_np <- vect("C:/Users/carme/NR6950/Final Project Part 2/Benoue/Benoue_NP/Benoue.shp")
## Warning: [vect] Z coordinates ignored
benoue <- project(benoue_np, crs(red_jan))
#cropping to Benoue NP polygon
red_list_cropped_masked <- lapply(red_list, function(r) {
r2 <- crop(r, benoue_33n)
mask(r2, benoue_33n)
})
sapply(red_list_cropped_masked, function(x) c(ext(x), ncol(x), nrow(x)))
## [,1]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,2]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,3]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,4]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,5]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,6]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
## [,7]
## [1,] <S4 class 'SpatExtent' [package "terra"] with 1 slot>
## [2,] 2216
## [3,] 2992
#Stacking
red_stack <- rast(red_list_cropped_masked)
names(red_stack) <- c("Jan","Mar","May","Jul","Sep","Oct","Nov")
#Checking projections of nir_stack and red_stack:
#crs(nir_stack)
#crs(red_stack)
#Calculating NDVI:
ndvi_stack <- (nir_stack - red_stack) / (nir_stack + red_stack)
## |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|========================================= |---------|---------|---------|---------|=========================================
#names(ndvi_stack) <- paste0("NDVI_", names(nir_stack))
plot(ndvi_stack)
#Writing each month to disk as GeoTIFF:
out_dir <- "C:/Users/carme/NR6950/Final Project Part 3/outputs"
for (i in 1:nlyr(ndvi_stack)) {
out_path <- file.path(out_dir, paste0("NDVI_", names(ndvi_stack)[i], ".tif"))
message("Writing: ", out_path)
writeRaster(ndvi_stack[[i]], out_path, overwrite = TRUE)
}
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Jan.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Mar.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_May.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Jul.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Sep.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Oct.tif
## Writing: C:/Users/carme/NR6950/Final Project Part 3/outputs/NDVI_Nov.tif
#Unsupervised Classification
library(terra)
k <- 10 # number of clusters
class_list <- list()
set.seed(99)
for (i in 1:nlyr(ndvi_stack)) {
cat("Processing layer:", names(ndvi_stack)[i], "\n")
# extract values for this single layer
vals <- values(ndvi_stack[[i]])
# find valid pixels
good <- !is.na(vals)
# run k-means ONLY on valid pixels
km <- kmeans(vals[good], centers = k, iter.max = 500, nstart = 5, algorithm = "Lloyd")
# create empty vector of all NA
out_vals <- rep(NA, length(vals))
# assign clusters to non-NA pixels
out_vals[good] <- km$cluster
# create classified raster for this layer
classified <- ndvi_stack[[i]]
classified[] <- out_vals
# store it
class_list[[i]] <- classified
}
## Processing layer: Jan
## Processing layer: Mar
## Processing layer: May
## Processing layer: Jul
## Processing layer: Sep
## Processing layer: Oct
## Processing layer: Nov
# convert list to SpatRaster stack
classified_stack <- rast(class_list)
names(classified_stack) <- paste0("Class_", names(ndvi_stack))
# Install if needed
#install.packages("scico")
# Create a color vector using the "navia" palette from scico
cols <- scico(10, palette = "navia") # n = number of classes or colors you want
# Plot your classified raster using these colors
plot(classified_stack, col = cols)
#Testing out leaflet code...
# Make sure classified_stack has layer names
#names(classified_stack) <- paste0("NDVI_", names(classified_stack))
# Create a color palette for the NDVI classes
# Assuming 10 classes (k = 10)
classes <- sort(unique(values(classified_stack)))
ndvi_pal <- colorFactor(
palette = scico::scico(length(classes), palette = "navia"),
domain = classes,
na.color = "transparent"
)
#ndvi_pal <- colorFactor(palette = scico(10, palette = "navia"), domain = unique(values(classified_stack), )
month_pal <- colorFactor(rainbow(12), points$month)
# Start Leaflet map
m <- leaflet() %>%
addProviderTiles("Esri.WorldImagery") %>%
# Add points
addCircleMarkers(
data = points,
lng = ~long, lat = ~lat,
color = ~month_pal(month),
radius = 4,
fillOpacity = 0.8,
stroke = FALSE,
popup = ~paste("Date:", timestamp),
group = "Points"
) %>%
addLegend(
"bottomright",
pal = month_pal,
values = points$month,
title = "Month"
)
# Add each NDVI layer ONCE
for(i in 1:nlyr(classified_stack)) {
m <- m %>%
addRasterImage(
classified_stack[[i]],
colors = ndvi_pal,
opacity = 1,
group = names(classified_stack)[i]
)
}
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
# Add a simple NDVI legend (one for all layers)
m <- m %>%
addLegend(
position = "bottomleft",
pal = ndvi_pal,
values = 1:10, # your k-means classes
title = "NDVI Classes",
opacity = 1,
labFormat = labelFormat(prefix = "Class ")
)
# Add layer control + hide unwanted layers
m <- m %>%
addLayersControl(
overlayGroups = c("Points", names(classified_stack)),
options = layersControlOptions(collapsed = FALSE)
) %>%
hideGroup(c("Class_Mar", "Class_May",
"Class_Jul", "Class_Sep",
"Class_Oct", "Class_Nov"))
m