ggplot()+geom_sf(data=esp0c, color="black", fill=NA)+ggtitle("The Largest Polygon associated with Spain")
Steps to obtain the figure
Here is how we select the largest polygon associated with Spain
#Step1: Loading librarieslibrary(sf)library(ggplot2)#Step2: Turn "s2" off to enable polygon simplificationsf_use_s2(FALSE)#Step3: Reading Spain bordersesp0 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet =TRUE)esp0<-st_simplify(esp0, dTolerance=0.01)#Note: Spain is made out of the polygons#We will focus on the largest polygons#that are traditionally association with Spain#Step4: Turning multipolygons into polygonsesp0b =st_cast(esp0,"POLYGON")#Step5: Calculating polygon areaesp0b$area<-st_area(esp0b)#Step6: Turning the units into square kmesp0b$area2<-units::set_units(esp0b$area, km^2)#Step7: Only choosing the largest polygonesp0c<-subset(esp0b, area2==max(esp0b$area2))#Step8: Making a box object around the polygonbbox <-st_bbox(esp0c)#Step9: Extracting min and max latitudes and longitudesmin_lat_y <- bbox["ymin"]max_lat_y <- bbox["ymax"]min_lon_x <- bbox["xmin"]max_lon_x <- bbox["xmax"]
ggplot()+geom_sf(data=esp0c, color="black", fill=NA)+ggtitle("The Largest Polygon associated with Spain")
Steps to obtain the figure
We will now process the fires data
#Step1: Reading wildfiresfires <-read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')#Step2: Turning the relevant time variable into a date variable recognized by Rfires$acq_date<-as.Date(fires$acq_date)#Step3: Turning the CSV file into an SF objectfires_sf<-st_as_sf(fires, coords =c("longitude", "latitude"))#Step4: Associating the crs information for the Spain shapefile to pointsfires_sf <-st_set_crs(fires_sf, st_crs(esp0)) # 4326 is the EPSG code for WGS84#Step5: Only selecting fires of type 0#See the variable definitions:#0 = presumed vegetation fire#1 = active volcano#2 = other static land source#3 = offshore detection (includes all detections over water)fires_sf<-subset(fires_sf, type==0)#Step6: Extract month informationfires_sf$month <-format(fires_sf$acq_date, "%m")#Step7: Reading only Augustfires_aug<-subset(fires_sf, acq_date>"2021-07-31"& acq_date<"2021-09-01")
#Defining an error for better map extenterror<-0.05ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_aug, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_viridis_c(option ="A", name ="bright_ti4")+scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +#We still need to define the extent: fires also happen outside the large polygon#associated with Spaincoord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Wild Fire Intensity in Spain, August '21") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))
Making the map
We can finally make the map:
Entire Year
We can also produce the following map
Entire Year
This is done using facet_wrap
ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_viridis_c(option ="A", name ="bright_ti4")+scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Wild Fire Intensity in Spain, August '21") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))+facet_wrap(~month, ncol =4) # Adjust the number of columns as needed
Entire Year
This is done using facet_wrap
ggplot() +geom_sf(data = esp0, color ="black", fill="black") +geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha =0.7) +scale_color_viridis_c(option ="A", name ="bright_ti4")+scale_size_continuous(range =c(0.5, 2), name ="bright_ti4") +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Wild Fire Intensity in Spain, August '21") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top"))+facet_wrap(~month, ncol =4) # Adjust the number of columns as needed
Entire Year
Or we can produce an animated map:
Entire Year
Or we can produce an animated map:
library("gganimate")#Extracting coordinatesfires_sf$lat_x<-st_coordinates(fires_sf)[,"X"]fires_sf$lon_y<-st_coordinates(fires_sf)[,"Y"]#Dropping the geometryfires<-st_drop_geometry(fires_sf)#Mappingfig<-ggplot() +geom_sf(data = esp0, color ="black", fill="black")+#We need to use geom_point NOT geom_sf because of library compatibility problemsgeom_point(data = fires, aes(x= lat_x, y = lon_y, color = bright_ti4, size = bright_ti4), alpha =0.9)+#Note that animated geom_sf with points does not work because of some library incompatibility#geom_sf(data = fires_sf, aes(color = bright_ti4, size = bright_ti4), alpha = 0.7) +scale_color_viridis_c(option ="A", name ="bright_ti4")+scale_size_continuous(range =c(0.1, 4), name ="bright_ti4") +guides(color =guide_legend(title ="bright_ti4"),size =guide_legend(title ="bright_ti4"))+coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error))+#animate argumentstransition_time(acq_date) +ggtitle('Date: {frame_time}')#defining animation options: speed of transition and resolutionanimate(fig, fps=2, res =60)
Spatial Join
spatial join is one of the most common operations in spatial analysis.
In a spatial join we are “attaching” attributes from one layer to another, based on their spatial relations.
For example, we will match the fires to Spanish provinces.
Spatial Join: Count by Polygon
Spatial Join: Count by Polygon
Place the folder with the shape files in the relevant week folder
#Step1: Reading wildfiresfires <-read.csv('./data/fires_2021/viirs-snpp_2021_Spain.csv')#Step2: Turning the acq_date variable into a date variablefires$acq_date<-as.Date(fires$acq_date)#Step3: Make sf object from pointsfires<-st_as_sf(fires, coords =c("longitude", "latitude"))fires <-st_set_crs(fires, st_crs(esp0)) # 4326 is the EPSG code for WGS84#Step4: Extract month and year informationfires$month <-format(fires$acq_date, "%m")fires$year <-format(fires$acq_date, "%Y")#Step5: Simplifying our datafires<-subset(fires, type==0)#Step6: Reading only Augustfires_aug<-subset(fires, acq_date>"2021-07-31"& acq_date<"2021-09-01")
Spatial Join: Count by Polygon
We will try to count the number of fires by province
Original Data Comparison
And this is how it differs from the original data: note the many overlapping points
Jittering the points allows us to see better the ones that are clustered
Show the code - Low Jitter
error<-0.05ggplot() +geom_sf(data = esp2b) +#low jittergeom_sf(data =st_jitter(fires_augb, 0.05), size=0.05, color ="red", alpha =0.7) +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Wild Fire Intensity in Spain, August '21") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top") )
Show the code - High Jitter
error<-0.05ggplot() +geom_sf(data = esp2b) +#high jittergeom_sf(data =st_jitter(fires_augb, 0.5), size=0.05, color ="red", alpha =0.7) +coord_sf(xlim =c(min_lon_x - error, max_lon_x + error), ylim =c(min_lat_y - error, max_lat_y + error)) +ggtitle("Wild Fire Intensity in Spain, August '21") +guides(color =guide_legend(title.position ="top"),size =guide_legend(title.position ="top") )
Original Data Comparison
And this is how the original data compares to the counts
Spatial Join: Count by Polygon - Step 1
And this is how we produce this graph:
Step 1:
library(dplyr)library(gridExtra)#Step1: Selecting only one variableesp2b<-subset(esp2, select =c("NAME_2"))#Step2: Selecting only bright_ti5. #This represents I-5 Channel brightness temperature of the fire pixel measured in Kelvin.fires_augb<-subset(fires_aug, select =c("bright_ti5"))#Step3: Adding a variable fire_no = 1 to help count firesfires_augb$fire_no<-1#Step4: Performing the spatial joinxjoin<-st_join(fires_augb, esp2b)
The unit is fire (i.e. the invidual event). The unit is NOT the district.
This is why Avila appears over 1,000 times
Thus, we need to count how many fires we have per district.
Spatial Join: Count by Polygon - Step 1
#Step4: Dropping the geometry from the object (i.e. turning it into a regular dataframe)#We need to drop the geometry to aggregate by groupxjoin<-st_drop_geometry(xjoin)#Step5: Summing by polygonssum_per_province <- xjoin %>%group_by(NAME_2) %>%summarize(sum_fires =sum(fire_no, na.rm =TRUE))#Step6: Performing a left join to the shapefilesesp1c<-left_join(esp2b, sum_per_province, by =c("NAME_2"="NAME_2"))#Step7: Replacing NAs with 0esp1c$sum_fires[is.na(esp1c$sum_fires)]<-0#Step8: Calculating logs to eliminate extreme valuesesp1c$log_fires<-log(esp1c$sum_fires+1)
Spatial Join: Count by Polygon - Step 2
And this is how we produce this graph:
Step 2:
#Step9: Calculating SD and Meanfires_sd <-sd(esp1c$log_fires, na.rm=T)fires_mean <-mean(esp1c$log_fires, na.rm=T)#Step10: Creating different stepsstep1<-min(esp1c$log_fires, na.rm=T)step2<-fires_mean - fires_sdstep3<-fires_mean -0.5*fires_sdstep4<-fires_meanstep5<-fires_mean +0.5*fires_sdstep6<-fires_mean + fires_sdstep7<-max(esp1c$log_fires, na.rm=T)#Step11: Putting all the steps in a listclint <-c(step1, step2, step3, step4, step5, step6, step7)clint
#Step1: Reading Spain borderssf_use_s2(FALSE)esp0 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_0.shp", quiet =TRUE)esp0<-st_simplify(esp0, dTolerance=0.05)#Step2: Simplifying country polygonesp0b =st_cast(esp0,"POLYGON")#Step3: Calculating polygon areaesp0b$area<-st_area(esp0b)#Step4: Turning the units into square kmesp0b$area2<-units::set_units(esp0b$area, km^2)#Step5: Only choosing the largest polygonesp0c<-subset(esp0b, area2==max(esp0b$area2))#Step6: Making a box object around the polygonbbox <-st_bbox(esp0c)#Step7: Extracting min and max latitudes and longitudesmin_lat_y <- bbox["ymin"]max_lat_y <- bbox["ymax"]min_lon_x <- bbox["xmin"]max_lon_x <- bbox["xmax"]
```{r , fig.width = 6.5, fig.height=2.4}library(gridExtra)#Arranging the Map and Histogram together#in grid.arrange there are 4 slots in a row. you need to indicate which slot will be taken#by fig1 vs. fig2. In our case, fig1 takes the first 3 slots#fig2 takes the 4th slot.pic<-grid.arrange(p1,p2, layout_matrix=rbind(c(1,1,1,2)))```
Spatial Join: Count by Polygon - Step 6
And this is how we produce this graph:
Step 6:
Spatial Join: Count by Polygon - Step 6
To save the map that you made still look good, you should use the fillowing dimensions:
#Step1: Reading Spain level-2 bordersesp2 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_2.shp", quiet =TRUE)#Step2: Simplifying polygonsesp2<-st_simplify(esp2, dTolerance=0.01)#Step3: Selecting only one variable of interestesp2<-subset(esp2, select=NAME_2)#Step4: Reading Spain level-4 borders - settlementsesp4 <-st_read(dsn="./data/gadm41_ESP_shp/gadm41_ESP_4.shp", quiet =TRUE)#Step5: Simplifying polygonsesp4<-st_simplify(esp4, dTolerance=0.01)#Step6: Selecting only one variable of interestesp4b<-subset(esp4, select=NAME_4)#Step7: Creating a unique ID (e.g. some settlements have the same name in multiple counties)esp4b$unique_id<-as.numeric(rownames(esp4b))#Step8: Calculating centroidesp4_pt<-st_centroid(esp4b)
We get the points esp4_pt (settlement centroids) to acquire the attributes of esp2 (county polygons)
#Step9: Perfoming the spatial joinpt_joined<-st_join(esp4_pt, esp2)#Step10: Dropping the geometry from the object (i.e. turning it into a regular dataframe)pt_joined_nogeom<-st_drop_geometry(pt_joined)
And this is the result:
#Step11: Inspecting the resulthead(pt_joined_nogeom, 10)
