library(sf)
library(ggplot2)
#Step1: Reading Spain borders
esp0 <- st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet = TRUE)
#Step2: Reading wildfires
fires <- read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')
#Step3: Turning the acq_date variable into a date variable
fires$acq_date<-as.Date(fires$acq_date)
#Step4: Make sf object from points
fires<-st_as_sf(fires, coords = c("longitude", "latitude"))
#Step5: Adding a CRS to the shapefile
fires <- st_set_crs(fires, st_crs(esp0))  # 4326 is the EPSG code for WGS84
#Step6: Extract month and year information
fires$month <- format(fires$acq_date, "%m")
fires$year <- format(fires$acq_date, "%Y")
#Step7: Simplifying our data
fires<-subset(fires, type==0)
#Step8: Reading only August
fires_aug<-subset(fires, acq_date>"2021-07-31" & acq_date<"2021-09-01")

#Step9: Simplify shapefiles
esp0b<-st_simplify(esp0, dTolerance=0.1)
#Step10: Simplifying country polygon
esp0 = st_cast(esp0b,"POLYGON")
#Step11: Calculating polygon area
esp0$area<-st_area(esp0)
#Step12: Turning the units into square km
esp0$area2<-units::set_units(esp0$area, km^2)
#Step13: Only choosing the largest polygon
esp0c<-subset(esp0, area2==max(esp0$area2))
#Step14: Making a box object around the polygon
bbox <- st_bbox(esp0c)
#Step15: Extracting min and max latitudes and longitudes
min_lat_y <- bbox["ymin"]
max_lat_y <- bbox["ymax"]
min_lon_x <- bbox["xmin"]
max_lon_x <- bbox["xmax"]

library(dplyr)
#Step16: Creating a grid around the country shapefile
r <- terra::rast(esp0, res=(0.5 * 0.5))
v <- terra::as.polygons(r, extent=FALSE)
#Step17: Turning the spatraster into an sf object
grid_sf<- sf::st_as_sf(v)
#Step18: Adding id
grid_sf$id<-rownames(grid_sf)
#Step19: Turning grids to points
grid_centr<-st_centroid(grid_sf)
#Step20: Performing the spatial join
grid_centr_joined<-st_join(grid_centr, esp0c)
#Step21: Dropping geometry
grid_centr_df1<-st_drop_geometry(grid_centr_joined)
#Step22: Dropping grids outside of the Spain polygons based on a variable
grid_centr_df2<-subset(grid_centr_df1, !is.na(COUNTRY))
#Step23: Selecting the id variable
grid_centr_df3<-subset(grid_centr_df2, select = c(id, COUNTRY))
#Step25: Joining grids with centroids
merged<-left_join(grid_sf, grid_centr_df3, by=c("id"="id"))
#Step26: Eliminating grids outside of the polygon
xsubset<-subset(merged, id %in% grid_centr_df3$id)

#Step27: Intersect points with grid
x <- st_intersects(fires_aug, xsubset, sparse = FALSE)
#Step28: Adding the points by grid
x <- colSums(x)
#Step29: Adding the sum the grid
xsubset$no_fires<-x
#Step30: Mapping the grid
fig2<-ggplot() +
  geom_sf(data=xsubset, aes(fill=no_fires)) + 
  theme_void() +
  scale_fill_viridis_b()+
  #geom_sf(data=esp0, fill=NA)+
  coord_sf(xlim = c(min_lon_x, max_lon_x), 
           ylim = c(min_lat_y, max_lat_y))

fig2

ggplot() +
  geom_sf(data=xsubset, aes(fill=log(no_fires+1))) + 
  theme_void() + scale_fill_viridis_b()+
  #geom_sf(data=esp0, fill=NA)+
  coord_sf(xlim = c(min_lon_x, max_lon_x), 
           ylim = c(min_lat_y, max_lat_y))