######"LavaFlow_mapper" version 1.0###### ###Tested in R 4.1.3 and RStudio 2022.02.01###### ###############CALL PACKAGES################ Be sure to install and call all the dependent packages before running library(tictoc) library(leaflet) library(leaflegend) library(sp) library(viridisLite) library(viridis) library(mapview) library(geosphere) library(htmlwidgets) #library(rgdal) #uncomment only if you want to plot a shapefile (WGS84) in the interactive thermal map. **For comparison purposes** tic() #start timer #################CLEAN VARIABLES AND TABS################## rm(list=ls()) #clean global environment graphics.off() #clean previous plots ####################SET WORKING DIRECTORY AND READ FIRMS DATASETS###################### setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) #set working directory automatically / this can be changed by R-users if needed SNPP<-read.csv(list.files(pattern = "SV(.*)csv$"),sep=",",dec=".",header=TRUE) #read Suomi data NOAA20<-read.csv(list.files(pattern = "J1V(.*)csv$"),sep=",",dec=".",header=TRUE) #read NOAA data #LF=readOGR(dsn=".",layer="Wolf_LF_22") #uncomment only if you have a shapefile with the area covered by the lavas; the name (here:Wolf_LF_22) must be replaced by your shapefile name in WGS84 ####################FILTER AND EXTRACT THE S-NPP DATA########################## Suomi<-subset(SNPP,track <= 0.5 & frp>=30) ## FILTERS!!! change the FRP filter according to the FRP_statistical results of your case study or use the default value of 35+-17 raw_times1=sprintf("%04d", unlist(Suomi["acq_time"]))#extract times and pad 0 if needed raw_dates1=sprintf(unlist(Suomi["acq_date"]))#extract dates datetime_raw1=Map(paste,raw_dates1,raw_times1)#merge dates and times datetime1=strptime(datetime_raw1,format = "%d/%m/%y %H%M")#dates and times as POSIX objects #datetime1=strptime(datetime_raw1,format = "%Y-%m-%d %H%M") #uncomment only if you have problems creating datetime1 FRP1=as.numeric(unlist(Suomi["frp"]))#extract FRP LATS1=as.numeric(unlist(Suomi["latitude"]))#extract latitudes LONGS1=as.numeric(unlist(Suomi["longitude"]))#extract longitudes DATE1=as.numeric(datetime1)/86400 #convert dates from seconds to days #########READ AND EXTRACT THE NOAA-20 DATA########### IMPORTANT: For past eruptions earlier than 2020 there is no NOAA-20 data, so we strongly recommend to duplicate part of the SNPP data to avoid errors. The duplicated data (3 rows) should surpass the track and FRP filters. NOAA<-subset(NOAA20,track <= 0.5 & frp>=30) ## FILTERS!!! change the FRP filter according to the FRP_statistical results of your case study or use the default value of 35+-17 raw_times2=sprintf("%04d", unlist(NOAA["acq_time"]))#extract times and pad 0 if needed raw_dates2=sprintf(unlist(NOAA["acq_date"]))#extract dates datetime_raw2=Map(paste,raw_dates2,raw_times2)#merge dates and times datetime2=strptime(datetime_raw2,format = "%d/%m/%y %H%M")#dates and times as POSIX objects #datetime2=strptime(datetime_raw2,format = "%Y-%m-%d %H%M") #uncomment only if you have problems creating datetime1 FRP2=as.numeric(unlist(NOAA["frp"]))#extract FRP LATS2=as.numeric(unlist(NOAA["latitude"]))#extract latitudes LONGS2=as.numeric(unlist(NOAA["longitude"]))#extract longitudes DATE2=as.numeric(datetime2)/86400 #convert dates from seconds to days ###############COLOR SCALE OF DATES################### seq<-c(min(min(DATE1),min(DATE2)):(max(max(DATE1),max(DATE2))+1)) #define min and max dates pal <- colorNumeric(palette = "Spectral", rev=TRUE,domain = seq) #use spectral palette for dates myLabelFormat = function(...,dates=FALSE){ if(dates){ function(type = "numeric", cuts){ as.Date(cuts, origin="1970-1-1")} } else {labelFormat(...) }} ## convert dates to numeric values to be used in the color-scale palette ###########PLOT THE THERMAL ANOMALIES ON AN ESRI BASEMAP############## if the interactive map doesn't appear in the viewer tab, select lines 56-64 and click on run m <- leaflet(Suomi) %>% addProviderTiles("Esri.WorldImagery") %>% #Add default basemaps (e.g. Esri.WorldImagery, OpenTopoMap, OpenStreetMap.Mapnik, for others see http://leaflet-extras.github.io/leaflet-providers/preview/index.html) addCircles(data=Suomi,lng=Suomi$longitude, lat=Suomi$latitude,weight = 3, radius=185, color = ~pal(DATE1), stroke = FALSE, fillOpacity = 0.7, popup = Suomi$acq_date) %>% #fillOpacity changes transparency addCircles(data=NOAA,lng=NOAA$longitude, lat=NOAA$latitude, weight = 3, radius=185, color = ~pal(DATE2), stroke = FALSE, fillOpacity = 0.7, popup = NOAA$acq_date) %>% #addPolygons(data = LF, fillOpacity = 0.2, color = "white",weight = 4)%>% #uncomment only if you want to compare with a map in shapefile format, LF=LavaFlow addLegend("topright", pal = pal, values = ~seq, title = "Date",opacity=1,labFormat = myLabelFormat(dates=TRUE))%>% #add color scale legend based on dates addScaleBar() #add scale bar m #print the map ##################SAVE AND EXPORT RESULTS################### if needed uncomment the line 66 #mapshot(m, file = "map.png") #save the map in png format; before using it, you must install the webshot package "webshot::install_phantomjs()". Otherwise you can export the map from the Viewer tab #saveWidget(m, file = "dynamic map.html") write.csv(Suomi,"filter_VIIRS_S-NPP.csv") #save the filtered S-NPP data write.csv(NOAA,"filter_VIIRS_NOAA-20.csv") #save the filtered NOAA-20 data ###############ESTIMATE LAVA FLOW DISTANCE IN KILOMETERS######### reference_long=-91.32241 #IMPORTANT!!! longitude coordinates must be replaced by those of your case study. If you do not have a reference point, keep the default, but the lenght estimations would be incorrect reference_lat=-0.01753 # IMPORTANT!!! latitude coordinates must be replaced by those of your case study. If you do not have a reference point, keep the default, but the lenght estimations would be incorrect distances1={}#initialize distances vector for S-NPP for (i in 1:length(LONGS1)){ distances1[i]=distm(c(reference_long, reference_lat), c(LONGS1[i], LATS1[i]))/1000#distance in kilometers for S-NPP } distances2={}#initialize distances vector for NOAA-20 for (i in 1:length(LONGS2)){ distances2[i]=distm(c(reference_long, reference_lat), c(LONGS2[i], LATS2[i]))/1000#distance in kilometers for NOAA-20 } ###############PLOT THE ESTIMATED LAVA FLOW DISTANCE######### maxx<-max(c(max(distances1),max(distances2))) #maximum distance between S-NPP and NOAA to define the vertical scale of the plot mindate<-as.Date(min(c(min(datetime1),min(datetime2)))) #minimum date between S-NPP and NOAA to define the horizontal scale of the plot maxdate<-as.Date(max(c(max(datetime1),max(datetime2)))) #maximum date between S-NPP and NOAA to define the horizontal scale of the plot plot(as.Date(datetime1),distances1, type='h',lwd=2,col="orange",xlab = "",xaxt="n",xlim=c(mindate,maxdate),ylim=c(0,maxx),ylab="Kilometers",main=sprintf("Estimated lava flow distances\n%s",paste(min(c(min(datetime1),min(datetime2))),max(c(max(datetime1),max(datetime2))),sep = ' to '))) #plot SNPP filtered data par(new=T) plot(as.Date(datetime2),distances2, type='h',lwd=2,col="purple3",xlab = "",xaxt="n",xlim=c(mindate,maxdate),ylim=c(0,maxx),ylab="") #plot NOAA filtered data axis(4, las=3) grid(nx = NA, ny=NULL, col = "lightgray", lty = "dotted",lwd = par("lwd")) #add grid axis.Date(1, at=seq(mindate, maxdate, by="week"), format="%b %d") #format of ticks legend("topleft", legend=c("S-NPP", "NOAA-20"),col=c("orange", "purple3"), lty=rep(1),lwd=2, cex=0.7, text.font=1, bty="o",box.lty = 1,bg = "white", ncol=2) #plot's legend df1 <- data.frame(Date_SNPP=c(datetime1),Distance_SNPP=c(distances1)) #create a dataframe of the distance per day based on S-NPP satellite df2<-data.frame(Date_NOAA=c(datetime2), Distance_NOAA=c(distances2)) #create a dataframe of the distance per day based on NOAA-20 satellite write.csv(df1,"distance_SNPP.csv") #save the distance based on S-NPP data write.csv(df2,"distance_NOAA.csv") #save the distance based on NOAA data toc() #end timer