Problem in applying 2D gaussian method to a raster

28 views Asked by At

I am applying the 2D Gaussian method to a raster image of LST but facing errors. the first plot is empty, opt has an infinite value, and no proper final output is shown.

gaussian2d <- function(x, y, mu_x, mu_y, sigma_x, sigma_y, amp)
{
  exp(-((x-mu_x)^2/(2*sigma_x^2) + (y-mu_y)^2/(2*sigma_y^2))) * amp
}

#define a function for the Getis Ord-Gi* statistic
getis_ord_gi <- function(x, y, z, w)
{
  n <- length(x)
  gi <- numeric(n)
  for(i in 1:n)
  { dist <- sqrt((x-x[i])^2 + (y-y[i])^2)
    neigh <- which(dist <= w)
    gi[i] <- (sum(z[neigh]) - mean(z) * length(neigh)) / (sd(z) * sqrt((n * length(neigh) - length(neigh)^2) / (n - 1)))
  }
  return(gi)
}
LST_Gi <- getis_ord_gi(coordinates(LST)[1], coordinates(LST)[2], values(LST), w=1000)
LST_Gi_raster <- rasterFromXYZ(cbind(coordinates(LST), LST_Gi))
plot(LST_Gi_raster, main="Getis-Ord-Gi* for LST")
Z_score <- LST_Gi/sd(LST_Gi)
p_value <- 2*pnorm(-abs(Z_score))
alpha <- 0.05
hot_spots <- which(p_value < alpha & Z_score > 0)
cold_spots <- which(p_value < alpha & Z_score < 0)
plot(LST, main = "Hot Spots and cold spots of UHI")
points(coordinates(LST)[hot_spots,], col="red", pch=16)
points(coordinates(LST)[cold_spots,], col="blue", pch=16)
legend("bottomright",legend=c("Hot Spots","Cold Spots"), col=c("red","blue"), pch=16)

x <- coordinates(LST)[hot_spots,1]
y <- coordinates(LST)[hot_spots,2]
z <- values(LST)[hot_spots]
mu_x0 <- mean(x)
mu_y0 <- mean(y)
sigma_x0 <- sd(x)/2
sigma_y0 <- sd(y)/2
amp0 <- max(z)

sse <- function(p)
{ mu_x <- p[1]
  mu_y <- p[2]
  sigma_x <- p[3]
  sigma_y <- p[4]
  amp0 <- p[5]
  fitted <- gaussian2d(x, y, mu_x, mu_y, sigma_x, sigma_y, amp)
  sum((z - fitted)^2)
}
opt <- optim(c(mu_x0, mu_y0, sigma_x0, sigma_y0, amp0), sse)
opt$par
persp(x, y, gaussian2d(x, y, opt$par[1], opt$par[2], opt$par[3], opt$par[4], opt$par[5]), theta = 30, phi = 20, col = "lightblue", main = "Fitted Gaussian surface")
threshold <- 0.5*opt$par[5]
uhi_extent <- which(gaussian2d(x, y, opt$par[1], opt$par[2], opt$par[3], opt$par[4], opt$par[5]) >= threshold)

plot(LST, main="Spatial extent of the urban heat island")
points(coordinates(LST)[hot_spots[uhi_extent],], col = "red", pch = 16) 
legend("bottomright", legend = "Spatial extent", col = "red", pch = 16)

ref_temp <- mean(values(LST)[cold_spots])
uhi_magnitude <- opt$par[5] - ref_temp

uhi_magnitude

I printed the parameters of optim function and learned these variables stored infinite values. I try to write those lines differently but the error is not resolved.

0

There are 0 answers