Package SPGWR:将GWR模型参数应用于更精细的空间比例

ttp71kqs  于 2022-12-20  发布在  其他
关注(0)|答案(1)|浏览(161)

我正在使用Rspgwr执行地理加权回归(GWR)。我想将模型参数应用于更精细的空间比例,但收到以下错误:Error in validObject(.Object): invalid class “SpatialPointsDataFrame” object: number of rows in data.frame and SpatialPoints don't match.
当我使用另一个名为GWmodel的GWR软件包时,我不会遇到这个问题。例如,使用GWmodel时,我会遇到:

# library(GWmodel)
# library(sp)
# library(raster)

ghs = raster("path/ghs.tif") # fine resolution raster
regpoints <- as(ghs, "SpatialPoints")

block.data = read.csv(file = "path/block.data.csv")

coordinates(block.data) <- c("x", "y")
proj4string(block.data) <- "EPSG:7767"

eq1 <- ntl ~ ghs
abw = bw.gwr(eq1, 
             data = block.data, 
             approach = "AIC", 
             kernel = "gaussian",
             adaptive = TRUE, 
             p = 2,
             parallel.method = "omp",
             parallel.arg = "omp")

ab_gwr = gwr.basic(eq1, 
                   data = block.data, 
                   regression.points = regpoints,
                   bw = abw, 
                   kernel = "gaussian", 
                   adaptive = TRUE, 
                   p = 2, 
                   F123.test = FALSE,
                   cv = FALSE,
                   parallel.method = "omp",
                   parallel.arg = "omp")

ab_gwr

sp <- ab_gwr$SDF
sf <- st_as_sf(sp)

# intercept
intercept = as.data.frame(sf$Intercept)
intercept = SpatialPointsDataFrame(data = intercept, coords = regpoints)
gridded(intercept) <- TRUE
intercept <- raster(intercept)
raster::crs(intercept) <- "EPSG:7767"

intercept = resample(intercept, ghs, method = "bilinear")
    
# slope
slope = as.data.frame(sf$ghs)
slope = SpatialPointsDataFrame(data = slope, coords = regpoints)
gridded(slope) <- TRUE
slope <- raster(slope)
raster::crs(slope) <- "EPSG:7767"

slope = resample(slope, ghs, method = "bilinear")

gwr_pred = intercept + slope * ghs

writeRaster(gwr_pred, 
            "path/gwr_pred.tif", 
            overwrite = TRUE)

如何使用spgwr包将GWR模型参数应用于更精细的空间比例?
下面是使用spgwr包的代码:

library(spgwr)
library(sf)
library(raster)
library(parallel)

ghs = raster("path/ghs.tif") # fine resolution raster
regpoints <- as(ghs, "SpatialPoints")

block.data = read.csv(file = "path/block.data.csv")

#create mararate df for the x & y coords
x = as.data.frame(block.data$x)
y = as.data.frame(block.data$y)

#convert the data to spatialPointsdf and then to spatialPixelsdf
coordinates(block.data) = c("x", "y")

# specify a model equation
eq1 <- ntl ~ ghs

# find optimal ADAPTIVE kernel bandwidth using cross validation
abw <- gwr.sel(eq1, 
               data = block.data, 
               adapt = TRUE, 
               gweight = gwr.Gauss)

# fit a gwr based on adaptive bandwidth
cl <- makeCluster(detectCores())
ab_gwr <- gwr(eq1, 
              data = block.data,
              adapt = abw, 
              gweight = gwr.Gauss,
              hatmatrix = TRUE,
              regpoints,
              predictions = TRUE,
              se.fit = TRUE,
              cl = cl)
stopCluster(cl)
#print the results of the model
ab_gwr

sp <- ab_gwr$SDF
sf <- st_as_sf(sp)

# intercept
intercept = as.data.frame(sf$Intercept)
intercept = SpatialPointsDataFrame(data = intercept, coords = regpoints)
gridded(intercept) <- TRUE
intercept <- raster(intercept)
raster::crs(intercept) <- "EPSG:7767"

intercept = resample(intercept, ghs, method = "bilinear")

# slope
slope = as.data.frame(sf$ghs)
slope = SpatialPointsDataFrame(data = slope, coords = regpoints)
gridded(slope) <- TRUE
slope <- raster(slope)
raster::crs(slope) <- "EPSG:7767"

slope = resample(slope, ghs, method = "bilinear")

gwr_pred = intercept + slope * ghs

writeRaster(gwr_pred, 
            "path/gwr_pred.tif", 
            overwrite = TRUE)

此外,如果我设置在

ab_gwr <- gwr(eq1, 
                  data = block.data,
                  adapt = abw, 
                  gweight = gwr.Gauss,
                  hatmatrix = TRUE,
                  fit.points = regpoints,
                  predictions = TRUE,
                  se.fit = TRUE,
                  cl = cl)

我收到此错误:Error in gwr(eq1, data = block.data, adapt = abw, gweight = gwr.Gauss,: No data slot in fit.points.
高分辨率光栅:ghs = raster(ncols=47, nrows=92, xmn=582216.388, xmx=603366.388, ymn=1005713.0202, ymx=1047113.0202, crs='+proj=lcc +lat_0=18.88015774 +lon_0=76.75 +lat_1=16.625 +lat_2=21.125 +x_0=1000000 +y_0=1000000 +datum=WGS84 +units=m +no_defs')
csv可以从here下载。

zfycwa2u

zfycwa2u1#

为了使用spgwr包将GWR的模型参数应用于更精细的空间比例:
1.以粗略比例计算GWR
1.使用参数fit.pointspredictionsfittedGWRobject再次应用步骤1。
代码:

library(spgwr)
library(sf)
library(raster)
library(parallel)

ghs = raster("path/ghs.tif") # fine res raster
regpoints <- as.data.frame(ghs[[1]], xy = TRUE)
regpoints = na.omit(regpoints)
colnames(regpoints)[3] = "the_name_of_the_fine_res_raster"
coordinates(regpoints) <- c("x", "y")
gridded(regpoints) <- TRUE

block.data = read.csv(file = "path/block.data.csv") # df containing the dependent and independent coarse variables

#convert the data to spatialPointsdf
coordinates(block.data) = c("x", "y")

# specify a model equation
eq1 <- ntl ~ ghs

# find optimal ADAPTIVE kernel bandwidth using cross validation
abw <- gwr.sel(eq1, 
               data = block.data, 
               adapt = TRUE, 
               gweight = gwr.Gauss)

# fit a gwr based on adaptive bandwidth
cl <- makeCluster(detectCores())
xx <- gwr(eq1, 
              data = block.data,
              adapt = abw, 
              gweight = gwr.Gauss,
              hatmatrix = TRUE,
              se.fit = TRUE,
              cl = cl)
stopCluster(cl)

# predict to a fine spatial scale
cl <- makeCluster(detectCores())
ab_gwr <- gwr(eq1, 
              data = block.data,
              adapt = abw, 
              gweight = gwr.Gauss,
              fit.points = regpoints,
              predictions = TRUE,
              se.fit = TRUE,
              fittedGWRobject = xx,
              cl = cl)
stopCluster(cl)

#print the results of the model
ab_gwr

sp <- ab_gwr$SDF
sf <- st_as_sf(sp)

# intercept
intercept = as.data.frame(sf$`(Intercept)`)
intercept = SpatialPointsDataFrame(data = intercept, coords = regpoints)
gridded(intercept) <- TRUE
intercept <- raster(intercept)
raster::crs(intercept) <- "EPSG:7767"

# slope
slope = as.data.frame(sf$ghs)
slope = SpatialPointsDataFrame(data = slope, coords = regpoints)
gridded(slope) <- TRUE
slope <- raster(slope)
raster::crs(slope) <- "EPSG:7767"

gwr_pred = intercept + slope * ghs

writeRaster(gwr_pred, 
            "path/gwr_pred.tif", 
            overwrite = TRUE)

相关问题