LPModel1.R

#Quick script for plotting variance measures of sky brightness for poster.
require(dplyr)
require(plyr)
require(ggplot2)
require(relaimpo)
require(matrixStats)
require(Hmisc)
require(corrplot)
require("PerformanceAnalytics")
require(ggmap)
require(maps)
require(mapview)
require(mapdata)
require(munsell)
require(leaflet)
require(devtools)
require(webshot)
require(stats)
require(segmented)
require(viridis)

wd <- "~/Desktop/Coastlight"
setwd(wd)

#Read in data generated by SQC batch analysis.
SQCData <- read.table("CoastlightSkyVarianceAllSites.csv", header=TRUE, sep=",",as.is=T,skip=0,fill=TRUE,check.names=FALSE, encoding = "UTF-8")
SQCData <- arrange(SQCData,`SQC File Name`)
#Rename verbose variables
colnames(SQCData)[which(names(SQCData) == "Scalar Illuminance")] <- "ScalarIlluminance"

#Calculate whole sky coefficient of variation in scalar illuminance
SkySectors <- (SQCData[grepl("D(.*?)R(.*?)Luminance", names(SQCData))])
SectorNames <- colnames(SkySectors)
SkySectors[SkySectors<0] <- 0
SkySectors$SDLuminance <- rowSds(as.matrix(SkySectors),na.rm=TRUE)
SkySectors$MeanLuminance <- rowMeans(SkySectors[SectorNames],na.rm=TRUE)
SkySectors$CoVLuminance <- SkySectors$SDLuminance/SkySectors$MeanLuminance
SkySectors$`SQC File Name` <- SQCData$`SQC File Name`
#Merge in these scalar illuminance statistics into the batch analyses data set.
SQCData <- left_join(SQCData,SkySectors[,c("SQC File Name","MeanLuminance","SDLuminance","CoVLuminance")],by=c("SQC File Name"))

#Read in 2012 VIIRS upwards radiance data per site.
VIIRS <- read.table("LPPointsUnder150.csv", header=TRUE, sep=",",as.is=T,skip=0,fill=TRUE,check.names=FALSE)
VIIRS[,'Latitude'] = round(VIIRS[,'Latitude'],5)
VIIRS[,'Longitude'] = round(VIIRS[,'Longitude'],5)
colnames(VIIRS)[which(names(VIIRS) == "Brightness (nW/Sr/cm^2)")] <- "VIIRSBrightness"

#Read in field data.
FieldData <- read.table("CoastlightFieldData.csv", header=TRUE, sep=",",as.is=T,skip=0,fill=TRUE,check.names=FALSE)
FieldData[,'Centroid latitude'] = round(FieldData[,'Centroid latitude'],5)
FieldData[,'Centroid longitude'] = round(FieldData[,'Centroid longitude'],5)
colnames(FieldData)[which(names(FieldData) == "SQM2015Atlas (mcd/m^2)")] <- "SQA"

#Merge in VIIRS data by pixel to field data.
FieldData <- left_join(FieldData,VIIRS[,c("Latitude","Longitude","VIIRSBrightness","OBJECTID")],by=c("Centroid latitude" = "Latitude","Centroid longitude" = "Longitude"))
#write.table(FieldData,"CoastligthFieldData.tsv",quote=FALSE,sep="\t",row.names = FALSE)

#Create unique site ID.  Each site contains 5 locations.
colnames(FieldData)[which(names(FieldData) == "OBJECTID")] <- "UniqueID"
FieldData$UniqueID <- as.factor(as.character(FieldData$UniqueID))

#Keep batch analysis data for images with zero horizon.
SQCDataZeroHorizon <- SQCData[!grepl("Horizon",SQCData$'SQC File Name'),]
SQCDataZeroHorizon <- arrange(SQCDataZeroHorizon,`SQC File Name`)

#Keep batch analysis data for images with edited horizon.
SQCDataEditedHorizon <- SQCData[grepl("Horizon",SQCData$'SQC File Name'),]
SQCDataEditedHorizon <- arrange(SQCDataEditedHorizon,`SQC File Name`)

#Zero-horizon merged data set
FieldSQCMergeZH <- left_join(FieldData,SQCDataZeroHorizon,by=c("SQCSiteName"="Location"))
#Edited-horizon merged data set
FieldSQCMergeEH <- left_join(FieldData,SQCDataEditedHorizon,by=c("SQCSiteName"="Location"))
#Determine the fraction of horizon radiance to full sky radiance.
FieldSQCMergeZH$HorizonLuminance <- (FieldSQCMergeZH$ScalarIlluminance - FieldSQCMergeEH$ScalarIlluminance)
FieldSQCMergeEH$HorizonLuminance <- (FieldSQCMergeZH$ScalarIlluminance - FieldSQCMergeEH$ScalarIlluminance)

#Subset horizon-edited data for analysis of spatial variability of scalar illuminance.
SpatialEH <- FieldSQCMergeEH
SpatialEH <- SpatialEH[,c("UniqueID","ScalarIlluminance")]
SpatialEH <- arrange(SpatialEH,UniqueID)
#Determine the mean of scalar illuminance by UniqueID
tmp <- as.data.frame(aggregate(SpatialEH$ScalarIlluminance,by=list(SpatialEH$UniqueID),FUN=mean))
colnames(tmp) <- c("UniqueID","ScalarIlluminanceSiteMean")
SpatialEH <- merge(SpatialEH,tmp,by=c("UniqueID"))
#Determine the coefficient of variation of scalar illuminance by UniqueID
tmp <- as.data.frame(aggregate(SpatialEH$ScalarIlluminance,by=list(SpatialEH$UniqueID),FUN=function(ScalarIlluminance){sd(ScalarIlluminance)/mean(ScalarIlluminance)}))
colnames(tmp) <- c("UniqueID","ScalarIlluminanceSiteCoV")
SpatialEH <- merge(SpatialEH,tmp,by=c("UniqueID"))
#Merge into larger data set.
FieldSQCMergeEH <- left_join(FieldSQCMergeEH,SpatialEH,by=c("UniqueID"))

#Subset zero horizon data for analysis of spatial variability of scalar illuminance.
SpatialZH <- FieldSQCMergeZH
SpatialZH <- SpatialZH[,c("UniqueID","ScalarIlluminance")]
SpatialZH <- arrange(SpatialZH,UniqueID)
#Determine the mean of scalar illuminance by UniqueID
tmp <- as.data.frame(aggregate(SpatialZH$ScalarIlluminance,by=list(SpatialZH$UniqueID),FUN=mean))
colnames(tmp) <- c("UniqueID","ScalarIlluminanceSiteMean")
SpatialZH <- merge(SpatialZH,tmp,by=c("UniqueID"))
#Determine the coefficient of variation of scalar illuminance by UniqueID
tmp <- as.data.frame(aggregate(SpatialZH$ScalarIlluminance,by=list(SpatialZH$UniqueID),FUN=function(ScalarIlluminance){sd(ScalarIlluminance)/mean(ScalarIlluminance)}))
colnames(tmp) <- c("UniqueID","ScalarIlluminanceSiteCoV")
SpatialZH <- merge(SpatialZH,tmp,by=c("UniqueID"))
#Merge into larger data set.
FieldSQCMergeZH <- left_join(FieldSQCMergeZH,SpatialZH,by=c("UniqueID"))

#Full merged data set
FieldSQCMergeEH$TypeHorizon <- "EditedHorizon"
FieldSQCMergeZH$TypeHorizon <- "ZeroHorizon"
FieldSQCMergeEH$`ScalarIlluminance.y` <- NULL
FieldSQCMergeZH$`ScalarIlluminance.y` <- NULL
colnames(FieldSQCMergeEH)[which(names(FieldSQCMergeEH) == "ScalarIlluminance.x")] <- "ScalarIlluminance"
colnames(FieldSQCMergeZH)[which(names(FieldSQCMergeZH) == "ScalarIlluminance.x")] <- "ScalarIlluminance"
#Use the same %horizon values for both the full hemispheric and edited horizon images.
FieldSQCMergeZH$Horizon <- FieldSQCMergeEH$Horizon
FieldSQCMerge <- rbind(FieldSQCMergeZH,FieldSQCMergeEH)
#write.table(FieldSQCMerge,"CoastlightProjectFinal.csv",quote=FALSE,sep=",",row.names = FALSE)

#World Atlas of the Artificial Night Sky Brightness = WAANSB
#Plot data for zero-horizon images
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=VIIRSBrightness,y=log10(ScalarIlluminance),color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=ScalarIlluminance))
LPPlotZH+xlab("VIIRS Upwards Radiance (nW/Sr/cm^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Cloud\ncover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=VIIRSBrightness,y=log10(ScalarIlluminance),color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=ScalarIlluminance))
LPPlotZH+xlab("VIIRS Upwards Radiance (nW/Sr/cm^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=ScalarIlluminance))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Cloud\ncover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=ScalarIlluminance))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx)")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=log10(ScalarIlluminance),y=CoVLuminance,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("Log(Scalar Illuminance (mlx))")+ylab("Coefficient of Variation\non Luminance")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=log10(ScalarIlluminance),y=CoVLuminance,color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("Log(Scalar Illuminance (mlx))")+ylab("Coefficient of Variation\non Luminance")+scale_color_gradientn("% Cloud\ncover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=VIIRSBrightness,y=ScalarIlluminanceSiteCoV,color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("VIIRS Upwards Radiance (nW/Sr/cm^2)")+ylab("Coefficient of Variation on\nScalar Illuminance within sites")+scale_color_gradientn("% Cloud\ncover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=VIIRSBrightness,y=ScalarIlluminanceSiteCoV,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("VIIRS Upwards Radiance (nW/Sr/cm^2)")+ylab("Coefficient of Variation on\nScalar Illuminance within sites")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=ScalarIlluminanceSiteCoV,color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Coefficient of Variation on\nScalar Illuminance within sites")+scale_color_gradientn("% Cloud\ncover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=ScalarIlluminanceSiteCoV,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=CoVLuminance))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Coefficient of Variation on\nScalar Illuminance within sites")+scale_color_gradientn("% Horizon",colours = plasma(10))
#

#Fit a linear model to predict the log of the scalar illuminance from field data for zero-horizon images.
LPModelZH <- lm(log10(ScalarIlluminance)~SQA+VIIRSBrightness,data=FieldSQCMergeZH)
summary(LPModelZH)
summary(aov(LPModelZH))
printCoefmat(coef(summary(step(LPModelZH))))
calc.relimp(LPModelZH)
#Initialize a simple linear model
LPLinearModelZH <- lm(ScalarIlluminance ~ SQA, data=FieldSQCMergeZH)

#Compare linear model residuals to cloud cover and horizon cover.
SQAIntercept <- as.numeric(lm(log10(FieldSQCMergeZH$ScalarIlluminance) ~ FieldSQCMergeZH$SQA)$coefficients[1])
SQASlope <- as.numeric(lm(log10(FieldSQCMergeZH$ScalarIlluminance) ~ FieldSQCMergeZH$SQA)$coefficients[2])
cor.test(FieldSQCMergeZH$Clouds,(SQAIntercept+SQASlope*FieldSQCMergeZH$SQA)-log10(FieldSQCMergeZH$ScalarIlluminance))
cor.test(FieldSQCMergeZH$Horizon,(SQAIntercept+SQASlope*FieldSQCMergeZH$SQA)-log10(FieldSQCMergeZH$ScalarIlluminance))

#Fit an asymptotic model to predict the log of the scalar illuminance from field data for zero-horizon images.
#SSasymp(input, Asym, R0, lrc)
#Asym+(R0-Asym)*exp(-exp(lrc)*input)
#input	a numeric vector of values at which to evaluate the model.
#Asym	a numeric parameter representing the horizontal asymptote on the right side (very large values of input).
#R0	a numeric parameter representing the response when input is zero.
#lrc	a numeric parameter representing the natural logarithm of the rate constant.
LPAssymptoticModelZH <- nls(log10(ScalarIlluminance)~SSasymp(SQA,Asym,R0,lrc),data=FieldSQCMergeZH)
summary(LPAssymptoticModelZH)
asymptote <- summary(LPAssymptoticModelZH)[10]$coefficients[1]
r_0 <- summary(LPAssymptoticModelZH)[10]$coefficients[2]
rate_constant <- summary(LPAssymptoticModelZH)[10]$coefficients[3]
#Calculate the best-fit line assuming an asymptotic model.
FieldSQCMergeZH$AsymptoticFit <-  asymptote+(r_0-asymptote)*exp(-exp(rate_constant)*FieldSQCMergeZH$SQA)
#Residuals of asymptotic fit versus cloud cover and horizon cover.
cor.test(FieldSQCMergeZH$Clouds,log10(FieldSQCMergeZH$ScalarIlluminance)-FieldSQCMergeZH$AsymptoticFit)
cor.test(FieldSQCMergeZH$Horizon,log10(FieldSQCMergeZH$ScalarIlluminance)-FieldSQCMergeZH$AsymptoticFit)

#Plot the asymptotic fit against the actual scalar illuminance.
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=log10(ScalarIlluminance),y=AsymptoticFit,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=AsymptoticFit))
LPPlotZH+xlab("Measured Log(Scalar Illuminance (mlx))")+ylab("Modeled Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=log10(ScalarIlluminance),y=AsymptoticFit,color=Clouds))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=AsymptoticFit))
LPPlotZH+xlab("Measured Log(Scalar Illuminance (mlx))")+ylab("Modeled Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Cloud\nCover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Clouds))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) asymptote+(r_0-asymptote)*exp(-exp(rate_constant)*x), size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Cloud\nCover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Horizon))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) asymptote+(r_0-asymptote)*exp(-exp(rate_constant)*x), size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=`CCT (cosine Corrected)`))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) asymptote+(r_0-asymptote)*exp(-exp(rate_constant)*x), size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("CCT (K)",colours = plasma(10))
#
#Plot differences between actual scalar illuminance and modeled scalar illuminance versus other factors.
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=Horizon,y=log10(ScalarIlluminance)-AsymptoticFit))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=log10(ScalarIlluminance)-AsymptoticFit))
LPPlotZH+xlab("% Horizon")+ylab("Measured-Modeled\nLog(Scalar Illuminance (mlx))")
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=`CCT (cosine Corrected)`,y=log10(ScalarIlluminance)-AsymptoticFit))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=log10(ScalarIlluminance)-AsymptoticFit))
LPPlotZH+xlab("CCT(K)")+ylab("Measured-Modeled\nLog(Scalar Illuminance (mlx))")
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=`CCT (cosine Corrected)`,y=log10(ScalarIlluminance)-AsymptoticFit,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=log10(ScalarIlluminance)-AsymptoticFit))
LPPlotZH+xlab("CCT(K)")+ylab("Measured-Modeled\nLog(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance)-AsymptoticFit,color=`CCT (cosine Corrected)`))+geom_point()+theme(text = element_text(size=25))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Measured-Modeled\nLog(Scalar Illuminance (mlx))")+scale_color_gradientn("CCT (K)",colours = plasma(10))


#Compare the asymptotic fit against the actual scalar illuminance for images above or below a certain threshold of horizon coverage.
FieldSQCMergeZHSubset <- subset(FieldSQCMergeZH,FieldSQCMergeZH$Horizon<14.9)
LPPlotZH <- ggplot(FieldSQCMergeZHSubset, aes(x=log10(ScalarIlluminance),y=AsymptoticFit,color=Horizon))+geom_point()+theme(text = element_text(size=25))+geom_smooth(method=glm, aes(fill=AsymptoticFit))
LPPlotZH+xlab("Measured Log(Scalar Illuminance (mlx))")+ylab("Modeled Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
cor.test(FieldSQCMergeZHSubset$AsymptoticFit,log10(FieldSQCMergeZHSubset$ScalarIlluminance))
#

#Try a logrithmic model of scalar illuminance versus WAANSB.
LPLogModelZH <- lm(log10(ScalarIlluminance)~log10(SQA),data=FieldSQCMergeZH)
FieldSQCMergeZH$LogFit <- log10(FieldSQCMergeZH$SQA*LPLogModelZH$coefficients[2])+LPLogModelZH$coefficients[1]
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Horizon))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) log10(x*LPLogModelZH$coefficients[2])+LPLogModelZH$coefficients[1], size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=Clouds))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) log10(x*LPLogModelZH$coefficients[2])+LPLogModelZH$coefficients[1], size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Cloud\nCover",colours = plasma(10))
#
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x=SQA,y=log10(ScalarIlluminance),color=`CCT (cosine Corrected)`))+geom_point()+theme(text = element_text(size=25))+stat_function(fun=function(x) log10(x*LPLogModelZH$coefficients[2])+LPLogModelZH$coefficients[1], size=2)
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Log(Scalar Illuminance (mlx))")+scale_color_gradientn("CCT (K)",colours = rev(plasma(10)))
#
#Plot differences between actual scalar illuminance and modeled scalar illuminance versus other factors.
LPPlotZH <- ggplot(subset(FieldSQCMergeZH,FieldSQCMergeZH$Horizon<=15 & FieldSQCMergeZH$`CCT (cosine Corrected)`<=3220), aes(x=SQA,y=log10(ScalarIlluminance)-LogFit,color=`CCT (cosine Corrected)`))+geom_point()+theme(text = element_text(size=25))
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Measured-LogFit\nLog(Scalar Illuminance (mlx))")+scale_color_gradientn("CCT (K)",colours = plasma(10))
#

#Fit a broken stick model of scalar illuminance to the Sky Quality Atlas.
#Initialize a simple linear model
LPLinearModelZH <- lm(log10(ScalarIlluminance) ~ SQA, data=FieldSQCMergeZH)
#Generate segmented model using an initial guess of breakpoints
LPSegmentedModelZH <- segmented(LPLinearModelZH, seg.Z = ~SQA,psi = list(SQA = c(2.5)))
summary(LPSegmentedModelZH)
#Get breakpoint using initial guess
LPSegmentedModelZH$psi
#Get slopes of segments
slope(LPSegmentedModelZH)
#Plot segmented model.
LPPlotZH <- ggplot(FieldSQCMergeZH, aes(x = SQA, y = log10(ScalarIlluminance), color = Horizon))+geom_point()+theme(text = element_text(size=25))+geom_line(data=FieldSQCMergeZH,aes(x=SQA,y=fitted(LPSegmentedModelZH)),color="black")
LPPlotZH+xlab("WAANSB (mcd/m^2)")+ylab("Modeled Log(Scalar Illuminance (mlx))")+scale_color_gradientn("% Horizon",colours = plasma(10))

##Correlations of variance light pollution variables.
#Each significance level is associated to a symbol : p-values(0, 0.001, 0.01, 0.05, 0.1, 1) <=> symbols(“***”, “**”, “*”, “.”, " “)
chart.Correlation(FieldSQCMergeZH[,c("Scalar Illuminance","SDLuminance","SQM2015Atlas (mcd/m^2)","Brightness (nW/Sr/cm^2)","Clouds","Horizon","HorizonLuminanceFraction")], histogram=FALSE, method="spearman")
chart.Correlation(FieldSQCMergeEH[,c("Scalar Illuminance","SDLuminance","SQM2015Atlas (mcd/m^2)","Brightness (nW/Sr/cm^2)","Clouds","Horizon","HorizonLuminanceFraction")], histogram=FALSE, method="spearman")

#To map various measures of coastal light pollution.
#MapCoordinates <- subset(FieldSQCMergeZH, FieldSQCMergeZH$SQA > 9.39)
MapCoordinates <- FieldSQCMergeZH
colnames(MapCoordinates)[which(names(MapCoordinates) == "Latitude")] <- "SQCLatitude"
colnames(MapCoordinates)[which(names(MapCoordinates) == "Longitude")] <- "SQCLongitude"
colnames(MapCoordinates)[which(names(MapCoordinates) == "Adjusted latitude")] <- "latitude"
colnames(MapCoordinates)[which(names(MapCoordinates) == "Adjusted longitude")] <- "longitude"
MapCoordinates <- MapCoordinates[!is.na(MapCoordinates$latitude) & !is.na(MapCoordinates$longitude),]
CalMap = leaflet(MapCoordinates) %>% 
  addTiles()
ColorScale <- colorNumeric(palette=plasma(10),domain=log10(FieldSQCMergeZH$ScalarIlluminance))
CalMap %>% addCircleMarkers(color = ~ColorScale(log10(ScalarIlluminance)), fill = TRUE,radius=2,fillOpacity = 0.1) %>% 
  addProviderTiles(providers$Esri.WorldTopoMap) %>%
  leaflet::addLegend(position="topright", pal=ColorScale,values=~log10(ScalarIlluminance),opacity=1,title="Log(Scalar Illuminance (mlx))")