code to extract MODIS fire data

MODIS_fire.R

library(rgeeExtra)
library(AnalystHelper)
library(EVERSpatDat)
library(rgee)
library(sf)
library(raster)
library(terra)
library(tmap)

# -------------------------------------------------------------------------
ee_Initialize(user="...",drive = TRUE); # insert uername
ee_users()
extra_Initialize(); #for rgeeExtra

# CRS
nad83.pro <- st_crs("EPSG:4269")
utm17 <- st_crs("EPSG:26917")
wgs84 <- st_crs("EPSG:4326")

tmap_mode("view")

# GIS ---------------------------------------------------------------------
# Fomr EVERSpatDat
data("EvPA")
data("nps_clipped")

bicy <- subset(nps_clipped,UNIT_CODE == "BICY")
bicy_enp <- subset(nps_clipped,UNIT_CODE%in%c("BICY","EVER"))

plot(st_geometry(bicy_enp))
AOI.buffer <- st_buffer(bicy_enp,1000)|>
  st_cast("POLYGON")|>
  st_union()|>
  st_as_sf()
AOI.buffer$region <- "BICY, ENP + 1km buff"
AOI.buffer$area_m2 <- st_area(AOI.buffer)|>as.numeric()
plot(st_geometry(AOI.buffer))

# ROI ---------------------------------------------------------------------
rlist <- st_bbox(st_transform(AOI.buffer,wgs84))
ROI <- c(rlist$xmin, rlist$ymin,
         rlist$xmax, rlist$ymin,
         rlist$xmax, rlist$ymax,
         rlist$xmin, rlist$ymax,
         rlist$xmin, rlist$ymin)
ROI.shp<-matrix(ROI, ncol = 2, byrow = TRUE) |>
  list() |>
  st_polygon() |>
  st_sfc() |>
  st_set_crs(4326)

ROI.shp2 <- ROI.shp|>st_transform(utm17)
plot(st_geometry(ROI.shp2))
plot(st_geometry(AOI.buffer),add=T)
st_centroid(ROI.shp)

ee_ROI <- ROI.shp |>
  st_transform(utm17)|>
  st_transform(wgs84)|>
  sf_as_ee()

ee_ever <- AOI.buffer|>
  st_transform(wgs84)|>
  sf_as_ee()


# MODIS data --------------------------------------------------------------
# Load MODIS burned area image collection and filter by date
# Using 2021 as an example case
dataset <- ee$ImageCollection('MODIS/061/MCD64A1')$
  filter(ee$Filter$date('2021-01-01', '2021-12-31'))$
  filterBounds(ee_ROI)

# Select the BurnDate band (also uncertainty)
burned_area <- dataset$select('BurnDate')
burned_area_img <- burned_area$mosaic(); # mosaic output to one image
burned_area_img <- burned_area_img$clip(ee_ROI); # clip to ROI

# Define visualization parameters
burned_area_vis <- list(
  min = 30.0,
  max = 341.0,
  palette = c('4e0400', '951003', 'c61503', 'ff1901')
)

# Center the map on southwest Florida and add the burned area layer
Map$setCenter(-81.13, 25.86, 9)
Map$addLayer(burned_area_img, burned_area_vis, 'Burned Area')

## Export burn area to google drive
task_img <- ee_image_to_drive(
  image = burned_area_img,                 # Your ee$Image
  description = 'burned_area_export',      # Task name
  folder = 'EarthEngine',                  # Drive folder name (must exist or will be created)
  fileNamePrefix = 'burned_area_2021',# Output file name (without extension)
  region = burned_area_img$geometry(),     # Export region
  scale = 500,                             # Spatial resolution in meters
  fileFormat = 'GeoTIFF'
)

# Start the export
task_img$start()

## To import from google drive
drive_find("burned_area_2021")
# Download to current directory
drive_download("EarthEngine/burned_area_2021_2025_07_15_15_16_44.tif",
               path="./GEE/burned_area_2021.tif",overwrite = TRUE)

# read from local folder
burned_area_2021 <- rast("./GEE/burned_area_2021.tif")|>
  project("EPSG:26917")
plot(burned_area_2021)

tm_shape(burned_area_2021)+tm_raster()

# Read in fire perimeters shapefile (from Malon et al)
# perim_shape <- sf::st_read(file.path(firehist_folder, "EVER_BICY_1978_2023_perim.shp"))
perim_shape <- sf::st_read("./data/EVER_BICY_1978_2023_perim.shp")

perim_shape_2021 <- subset(perim_shape,Year == 2021)|>
  st_transform(utm17)|>
  subset(!is.na(Fir_Typ))

bbox.lims <- st_bbox(AOI.buffer)
bks.vals <- seq(1,140,30)
pal <- viridis::viridis(length(bks.vals)-1,alpha=0.75)
plot(st_geometry(perim_shape_2021),ylim=bbox.lims[c(2,4)],xlim=bbox.lims[c(1,3)],border="red",col="pink")
plot(st_geometry(bicy_enp),add=T,border="dodgerblue1")
image(burned_area_2021,add=T,breaks=bks.vals,col = pal)