[R-sig-Geo] Find the optimal model using step wise linear regression
I have a sf object with 2 columns, named fclass and geometry. I am
interested in the fclass column. The fclass has several unique values. I
convert this sf object to spatraster, I resample it and I perform a linear
regression where I use another spatraster, called ntl, as my response
(lm(ntl ~ road)).
My goal is to find which unique values give the largest r-squared in the lm
model using stepwise linear regression whith backward elimination. My
initial try was:
library(pacman)
pacman::p_load(terra, sf, dplyr)
ntl <- rast("path/ntl.tif") # response
v <- st_read("path/road17.shp") # sf object
# baseline r2
vterra <- vect(v)
ref <- rast("path/pop.tif") # get ext and pixel size
r <- rast(v, res = res(ref), ext = ext(ref))
x <- rasterizeGeom(vterra, r, "length", "m")
x_res <- resample(x, ntl, "average")
s <- c(ntl, x_res)
names(s) <- c("ntl", "road")
m <- lm(ntl ~ road, s, na.action = na.exclude)
baseline <- sqrt(summary(m)$adj.r.squared)
orig_baseline <- sqrt(summary(m)$adj.r.squared)
classes <- unique(v$fclass)
inclasses <- unique(v$fclass)
i <- 1
while (i <= length(classes)) {
class <- classes[i]
print(paste0("current class ", class))
print(paste0("orig baseline: ", orig_baseline, " - baseline: ", baseline))
print(classes)
v_filtered <- v[v$fclass != class, ]
vterra <- vect(v_filtered)
r <- rast(v, res = res(ref), ext = ext(ref))
x <- rasterizeGeom(vterra, r, "length", "m")
x_res <- resample(x, ntl, "average")
s <- c(ntl, x_res)
names(s) <- c("ntl", "road")
m <- lm(ntl ~ road, s, na.action = na.exclude)
class_r2 <- sqrt(summary(m)$adj.r.squared)
if (class_r2 > baseline) {
classes <- classes[-i]
baseline <- class_r2
} else {
print("class_r2 is less than baseline")
print(paste0("class_r2: ", class_r2, " - baseline: ", baseline))
i <- i + 1
}
}
but it's not efficient (it takes ~ 5 minutes).
I have found the MASS package which has the lmStepAIC function but I am not
sure how to select the unique values from the fclass column to use as
predictiors.
head(v, 6)
Simple feature collection with 6 features and 1 field
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 598675.9 ymin: 7111459 xmax: 609432.8 ymax: 7118729
Projected CRS: WGS 84 / UTM zone 35S
fclass geometry
1 secondary MULTILINESTRING ((598675.9 ...
2 secondary MULTILINESTRING ((600641.7 ...
3 residential MULTILINESTRING ((601734.8 ...
4 residential MULTILINESTRING ((601163.9 ...
5 residential MULTILINESTRING ((601104.2 ...
6 motorway_link MULTILINESTRING ((609432.8 ...
The unique values in the fclass column:
unique(v$fclass)
[1] "secondary" "residential""motorway_link" "service"
"primary""unclassified" "motorway"
[8] "tertiary" "trunk" "primary_link" "footway"
"track" "secondary_link" "unknown"
[15] "living_street" "pedestrian" "path" "bridleway"
"steps" "trunk_link" "track_grade1"
[22] "track_grade3" "track_grade5" "cycleway" "track_grade4"
"tertiary_lin
The dataset can be downloaded from my GoogleDrive
<https://drive.google.com/drive/folders/1a0Khhej4koA0uhYTrQmqmMg02qfRWbjt?usp=drive_link>
.
Thank you.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Plot the location of a break point identified by the bfastlite function
I am using the bfastlite funtion from the bfast package for time-series
analysis of nighttime light (NTL) imagery. I have monthly data from April
2013 to December 2022 (129 monthly images). Using the code below, the
bfastlite found 2 breaks in the time-series data:
library(bfast)
library(terra)
wd <- "path/"
l <- list.files(paste0(wd), pattern = '.tif$',
all.files = TRUE, full.names = FALSE)
rr <- lapply(paste0(wd, l), rast)
standard <- rr[[1]]
rs <- list(standard)
for (i in 2:length(rr)) {
rr[[i]] <- terra::crop(rr[[i]], standard, extend = TRUE)
}
s <- rast(rr)
m <- terra::as.matrix(s)
m <- m[!rowSums(is.na(m)), ]
# convert the matrix to timeseries matrix
tm <- ts(data = m, start = c(2013, 4), end = c(2022, 12), frequency = 12,
class = "ts")
bf <- bfastlite(
tm,
formula = response ~ harmon,
order = 3,
breaks = "all",
lag = 17,
slag = 7,
na.action = na.omit,
stl = "none",
decomp = c("stl", "stlplus"),
sbins = 1,
h = .25,
level = 0.5,
type = c("OLS-MOSUM", "RE"),
plot = TRUE,
hpc = "foreach"
)
bf
plot(bf)
How can I plot the location (i.e., pixels) of the breakpoints? You can
download the dataset from my GoogleDrive
<https://drive.google.com/drive/folders/1UloTKYHa7I3Wi_2pt3n3xbZ6vFvVz3_Y?usp=sharing>.
The entire code runs in less than 10 seconds.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Error when running bfastlite using a monthly time-series matrix as input
I have 12 raster images which I stack them. Then I converted the rasterstack to a matrix, I remove the NA's from the matrix and I convert the matrix to monthly time-series matrix. Finally, I run the bfastlite function, like so: library(bfast) library(raster) wd <- "path" l <- list.files(paste0(wd), pattern='.tif$', all.files=TRUE, full.names=FALSE) r <- raster::stack(paste0(wd, l)) m <- raster::as.matrix(r) m <- m[!rowSums(is.na(m)), ] # convert the matrix to timeseries matrix tm <- ts(data = m, start = c(2019, 1), end = c(2019, 12), frequency = 12, class = "ts") class(tm) bf <- bfastlite( tm, formula = response ~ trend + harmon, order = 3, breaks = "LWZ", na.action = na.omit, stl = "trend", decomp = "stl", sbins = 1, h = 0.5, level = 0.5, type = "OLS-MOSUM" ) And the error: Error in stats::stl(x, s.window = "periodic") : series is not periodic or has less than two periods I have tried many variations in the parameters of the bfastlite function. Any ideas what that error means and how I can overcome it? The dataset can be downloaded from my GoogleDrive ( https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7). RStudio 2023.12.0+369, R 4.3.2, Windows 11. [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Preprocess daily Black Marble (NVP46A2) nighttime light product
I downloaded NASA's Black Marble daily product (VNP46A2) which is in h5 format (the data can be dowloaded from their website <https://ladsweb.modaps.eosdis.nasa.gov/search/order/1/VNP46A2--5000> (it's free you just need to create and account) or from my GDrive <https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7>). One needs to preprocess the data using the Scientific Data Sets (SDS) included in the h5 file. Based on the User Guide, these are the following parameters I need to account for: Table 4, page 14, Value of QF_Cloud_Mask in the VNP46A1/VJ146A1 product: Bit -- Flag description key -- Interpretation 0 -- Day/night -- 0 = Night 4-5 -- Cloud Mask Quality -- 11 = High 6-7 -- Cloud Detection Results & Confidence Indicator -- 00 = Confident Clear 8 -- Shadow Detected -- 0 = No 9 -- Cirrus Detection (IR) (BTM15 – BTM16) -- 0 = No Cloud 10 -- Snow/ Ice Surface -- 0 = No Snow/Ice Table 7, page 17, Values of the Mandatory_Quality_Flag in VNP46A2/VJ146A2 product: Value -- Retrieval quality -- Algorithm instance 00 -- High-quality -- Main algorithm (Persistent nighttime lights) Table 8, page 17, Values of the Snow_Flag in VNP46A2/VJ146A2 product: Flag description key -- Value -- Interpretation Snow/ Ice Surface -- 00 -- No Snow/Ice In the above tables I included the bit values I want to use to preprocess the NTL product, called Gap_Filled_DNB_BRDFCorrected_NTL. I am using R's terra package to preprocess the product. So far what I have managed to do is: library(terra) wd <- "path/" r <- rast(paste0(wd, "VNP46A2.A2018038.h28v07.001.2020333204506.h5")) crs(r) <- "epsg:4326" # dimensions 2400*(15/(60*60)) h = 28 v = 7 ext(r) = c(-180+h*10,-180+(h+1)*10, (8-v)*10,(8-v+1)*10) # up to this point the code works well # the tif images inside the h5 file (for the ifel function below) ntl <- r[[3]] # this is the Gap_Filled_DNB_BRDFCorrected_NTL latest_high_quality_retrieval <- r[[4]] mandatory_quality_flag <- r[[5]] qf_cloud_mask <- r[[6]] snow_mask <- r[[7]] # here is the issue!!! result <- ifel(r[[4]] > 0 & r[[5]] == 00 & r[[6]] == 1 & r[[7]] == 00, r[[3]], NA) # scale factor based on the User Guide table 6, page 16 result1 <- result * 0.1 writeRaster(result1, paste0(wd, "ntl.tif"), overwrite = TRUE) The writeRaster function returns an empty raster with null values. I can understand how to syntax the ifel for the Mandatory_Quality_Flag and Snow_Flag but I can't understand how to syntax the bits for the QF_Cloud_Mask. Could you help me syntax the ifel function properly using the bits from the tables? In a very abstract sense, the ifel statement should say: If snow_flag is 00 AND Mandatory_Quality_Flag is 00 AND the bit 0 from the QF_Cloud_Mask is 0 AND the bit 4-5 from the QF_Cloud_Mask is 11 AND the bit 6-7 from the QF_Cloud_Mask is 0 AND the bit 8 from the QF_Cloud_Mask is 0 AND the bit 9 from the QF_Cloud_Mask is 0 AND the bit 10 from the QF_Cloud_Mask is 0 THEN keep the values of the Gap_Filled_DNB_BRDFCorrected_NTL ELSE NA -- Tziokas Nikolaos Cartographer Tel:(+44)07561120302 LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130> [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Spatial structure in the random forest (regression) residuals
I have a response variable which is the nighttime light luminosity (called
*ntl* in the data.frame I am using). As you can see in the image attached
(link below), it has some bright spots (highlighted in red) which are areas
with high brightness (outliers in the data.frame). Also, I shared the
histogram of the response. It's clear that the distribution is right-skewed
(and possibly there is bimodality?).
*Purpose of my analysis*
My goal is to predict the *ntl* from the coarse spatial resolution to a
higher. This means that, I need to maintain these bright spots (the
outliers) in the predicted image. Because at a later stage, I will
downscale the residuals of the regression using area-to-point Kriging, I
need the residuals to be random (no spatial structure).
*Analysis*
Following this <https://juliasilge.com/blog/ikea-prices/>approach, the
resulting RF residuals show clearly a spatial structure (i.e., the
residuals do not show a random pattern), especially it seems that the model
cannot capture the areas with high luminosity. This can also be seen if I
plot the fitted vs residuals (please see the GDrive link below). Again, for
a "good" model I would like to see a symmetric scatter of points around the
horizontal like at zero, indicating random deviations of predictions from
the observed values, but from the plot above I see that this isn't the case.
My thought is that in the NTL image one can see where the bright areas are
and were not, while in the covariates there is no such "variation" in these
areas. This means that, the RFR model I use cannot capture this variability
(I'm just thinking out loud).
What are your thoughts on that? How can I improve the model in order to
eliminate the spatial structure in residuals while mainting the bright
areas in the predicted high resolution image? Here is the complete code:
wd <- "path/"
# this a projected CRS
provoliko <- "EPSG:24313"
# the data.frame of the variables
df <- read.csv(paste0(wd, 'block.data.csv'))
# here I extract the coordinates
crds <- df[, 1:2]
# here I remove the coordinates from the data.frame
df <- df[, 3:12]
# random forest regression
library(tidymodels)
set.seed(456)
# splitting the data.frame into training and test set using stratified
random sampling
ikea_split <- initial_split(df, strata = ntl, prop = 0.80)
ikea_train <- training(ikea_split)
ikea_test <- testing(ikea_split)
set.seed(567)
ikea_folds <- bootstraps(ikea_train, strata = ntl)
ikea_folds
library(usemodels)
eq1 <- ntl ~ .
use_ranger(eq1, data = ikea_train)
library(textrecipes)
ranger_recipe <-
recipe(formula = eq1, data = ikea_train)
ranger_spec <-
rand_forest(mtry = tune(), min_n = tune(), trees = 1001) %>%
set_mode("regression") %>%
set_engine("ranger")
ranger_workflow <-
workflow() %>%
add_recipe(ranger_recipe) %>%
add_model(ranger_spec)
set.seed(678)
doParallel::registerDoParallel()
ranger_tune <-
tune_grid(ranger_workflow,
resamples = ikea_folds,
grid = 10
)
show_best(ranger_tune, metric = "rsq")
autoplot(ranger_tune)
final_rf <- ranger_workflow %>%
finalize_workflow(select_best(ranger_tune, metric = "rsq"))
final_rf
ikea_fit <- last_fit(final_rf, ikea_split)
ikea_fit
collect_metrics(ikea_fit)
collect_predictions(ikea_fit) %>%
ggplot(aes(ntl, .pred)) +
geom_abline(lty = 2, color = "gray50") +
geom_point(alpha = 0.5, color = "midnightblue") +
coord_fixed()
predict(ikea_fit$.workflow[[1]], ikea_test[15, ])
ikea_all <- predict(ikea_fit$.workflow[[1]], df[, 2:11])
ikea_all <- as.data.frame(cbind(df$ntl, ikea_all))
ikea_all$rsds <- ikea_all$`df$ntl` - ikea_all$.pred
ikea_all <- as.data.frame(cbind(crds, ikea_all))
ikea_all <- ikea_all[,-3:-4]
# export the residuals as a raster
library(terra)
rsds <- rast(ikea_all, type = "xyz")
crs(rsds) <- provoliko
writeRaster(rsds, paste0(wd, "rsds.tif"), overwrite = TRUE)
library(vip)
imp_spec <- ranger_spec %>%
finalize_model(select_best(ranger_tune, metric = "rsq")) %>%
set_engine("ranger", importance = "permutation")
workflow() %>%
add_recipe(ranger_recipe) %>%
add_model(imp_spec) %>%
fit(ikea_train) %>%
extract_fit_parsnip() %>%
vip(aesthetics = list(alpha = 0.8, fill = "midnightblue"))
Because the csv I'm using has several thousands of rows, I can share it via
a link, from here
<https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7?usp=sharing>
(also in the link you can see the spatial structure of residuals and
histogram of the response variable). Just so you know, running and tuning a
model takes ~ 2-3 minutes on my laptop (8 gigs of RAM, 4 cores processor
(CPU: Intel(R) Core(TM) i5-4210U CPU @ 1.70GHz)). i am using R 4.3.1 and
RStudio 2023.09.0
[R-sig-Geo] NASA's Black Marble monthly data: Reprojection isn't accurate
I downloaded NASA's Black Marble monthly nighttime light NTL data, VNP46A3 <https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A3/>. In a previous question <https://gis.stackexchange.com/questions/466571/extent-not-found-on-nasas-black-marble-monthly-images-how-to-set-it/466574?noredirect=1#comment761916_466574> of mine the reprojection worked perfectly but now it seems that it doesn't. For example, I wanted to download NTL data for the city of Mumbai, India. After reprojecting the NTL (product 5 (All_Angles_Snow_Free) from the .h5) the result is attached here <https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7>. At the bottom if the image is a shp of Mumbai (downloaded from GADM) and the red circle in the top indicates where Mumbai is in the NTL image. Clearly something's not right. I downloaded the image from here <https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A3/> (LAADS-DAAC, Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center). The code I used to extract the NTL radiance image is: library(terra) wd <- "path/" r <- rast(paste0(wd, "VNP46A3.A2018091.h25v07.001.2021125122857.h5")) crs(r) <- "epsg:4326" 2400*(15/(60*60)) h = 25 v = 7 ext(r) = c(-180+h*10,-180+(h+1)*10, (v-2)*10,(v-1)*10) ntl <- r[[5]] writeRaster(ntl, paste0(wd, "ntl.tif"), overwrite = TRUE) Why the code worked perfectly in the previous question and now it doesn't? >From here <https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7> you can download the .h5 image if you don't want to use NASA's website. I am using R 4.3.1 and RStudio 2023.06.2+561. -- Tziokas Nikolaos Cartographer Tel:(+44)07561120302 LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130> [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] plot_tiltedmaps isn't showing anything
I am trying to plot 2 raster files but when I am using the *plot_tiltedmaps
*function the plot is empty. I am following the steps in the 2nd tutorial
<https://github.com/marcosci/layer>, called *More advanced example*. Any
ideas why is that?
Here is the code:
# install.packages("remotes")
remotes::install_github("marcosci/layer")
library(layer)
library(terra)
wd <- "path/"
pop = rast(paste0(wd, "pop.tif"))
agbh = rast(paste0(wd, "agbh.tif"))
pop <- tilt_map(pop, y_tilt = 3, x_shift = 25, y_shift = 50, parallel =
TRUE)
agbh <- tilt_map(agbh, y_tilt = 3, x_shift = 50, y_shift = 100, parallel =
TRUE)
map_list <- list(pop, agbh)
plot_tiltedmaps(map_list, palette = c("tofino", "rocket"), direction =
c(-1, 1))
You can download the data from here
<https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7?usp=drive_link>
.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Spatial random forest prediction: Error when predicting at unseen locations at finer spatial scale
I am using the package spatialRF in R for a spatial random forest
regression (SRFR) task. I have one response variable and 4 predictors and I
am performing SRFR at a coarse spatial scale. My goal is to take the model
parameters and apply them to a finer spatial resolution in order to predict
the response variable at the finer spatial scale.
When I run
p <- stats::predict(object = model.spatial, #name of the
spatialRF model
data = s, # data.frame
containing the predictors at the fine spatial scale (without NaN values)
type = "response")$predictions
I am getting this error: Error in predict.ranger.forest(forest, data,
predict.all, num.trees, type,: Error: One or more independent variables not
found in data.
I have checked the column names of s and my original data.frame (the one I
used to build the model at the coarse scale) and they are the same. How can
I use the model i created at the coarse scale to predict the response
variable at a finer spatial scale?
Here is the code:
library(spatialRF)
library(stats)
wd = "path/"
block.data = read.csv(paste0(wd, "block.data.csv")) # coarse resolution
#names of the response variable and the predictors
dependent.variable.name <- "ntl"
predictor.variable.names <- colnames(block.data)[4:7]
#coordinates of the cases
xy <- block.data[, c("x", "y")]
block.data$x <- NULL
block.data$y <- NULL
#distance matrix
distance.matrix <- as.matrix(dist(block.data))
min(distance.matrix)
max(distance.matrix)
#distance thresholds (same units as distance_matrix)
distance.thresholds <- c(0, 20, 50, 100, 200, 500)
#random seed for reproducibility
random.seed <- 456
#creating and registering the cluster
local.cluster <- parallel::makeCluster(
parallel::detectCores() - 1,
type = "PSOCK")
doParallel::registerDoParallel(cl = local.cluster)
# fitting a non-spatial Random Forest
model.non.spatial <- spatialRF::rf(
data = block.data,
dependent.variable.name = dependent.variable.name,
predictor.variable.names = predictor.variable.names,
distance.matrix = distance.matrix,
distance.thresholds = distance.thresholds,
xy = xy,
seed = random.seed,
verbose = FALSE)
# Fitting a spatial model with rf_spatial()
model.spatial <- spatialRF::rf_spatial(
model = model.non.spatial,
method = "mem.moran.sequential",
verbose = FALSE,
seed = random.seed)
#stopping the cluster
parallel::stopCluster(cl = local.cluster)
# prediction at a finer spatial scale
s = read.csv(paste0(wd, "s.csv")) # df containg the predictors at fine
scale
p <- stats::predict(object = model.spatial,
data = s,
type = "response")$predictions
I tried solutions like:
levels(s$lc) <- levels(block.data$lc)
in case I had missing land cover types in the lc column between the spatial
scales, but it didn't work.
>From here
you can download the two data.frames.
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Extract residuals from an XGBoost regression as a single raster
952514648, 6.28353786468506, 4.67424774169922,
4.56812858581543, 1.88153350353241, 4.31531000137329)), class =
"data.frame", row.names = c(NA,
-34L))
My raster layer:
r = new("RasterBrick", file = new(".RasterFile", name = "", datanotation =
"FLT4S",
byteorder = "little", nodatavalue = -Inf, NAchanged = FALSE,
nbands = 1L, bandorder = "BIL", offset = 0L, toptobottom = TRUE,
blockrows = 0L, blockcols = 0L, driver = "", open = FALSE),
data = new(".MultipleRasterData", values = structure(c(NA,
NA, NA, NA, 18.7969169616699, NA, NA, NA, NA, NA, 25.7222957611084,
23.4188251495361, 25.4322757720947, 16.4593601226807, 12.7868213653564,
NA, NA, NA, NA, 30.9337253570557, 29.865758895874, 30.4080600738525,
29.5479888916016, 24.3493347167969, NA, NA, NA, 35.2427635192871,
38.989933013916, 34.6536979675293, 29.4607238769531, 30.7469024658203,
34.3946380615234, NA, NA, NA, NA, 42.8660278320312, 34.7930717468262,
30.9516315460205, 32.20654296875, 39.999755859375, 46.6002235412598,
38.6480979919434, NA, NA, 60.5214920043945, 33.1799964904785,
31.8498134613037, 30.9209423065186, 32.2269744873047, 53.7062034606934,
45.5225944519043, 38.3570976257324, NA, 123.040382385254,
73.0528182983398, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
19.6407718658447, NA, NA, NA, NA, NA, 610.009216308594,
654.812622070312,
426.475830078125, 66.3839492797852, 10.6471328735352, NA,
NA, NA, NA, 443.848846435547, 602.677429199219, 488.478454589844,
387.470947265625, 58.2341117858887, NA, NA, NA, 413.888488769531,
315.057678222656, 354.082946777344, 602.827758789062, 463.518829345703,
296.713928222656, NA, NA, NA, NA, 923.920593261719, 434.436645507812,
799.562927246094, 404.709564208984, 265.043304443359, 366.697235107422,
399.851684570312, NA, NA, 952.2314453125, 870.356994628906,
673.406616210938, 493.521606445312, 273.841888427734, 371.428619384766,
383.057830810547, 320.986755371094, NA, 991.131225585938,
1148.87768554688, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
39.7242431640625, NA, NA, NA, NA, NA, 44.9583969116211,
41.4048385620117,
42.6056709289551, 40.0976028442383, 38.7490005493164, NA,
NA, NA, NA, 44.2747650146484, 43.5645370483398, 41.6180191040039,
40.3799781799316, 38.8664817810059, NA, NA, NA, 44.9089202880859,
44.414306640625, 44.560977935791, 43.1288986206055, 40.9315185546875,
38.8918418884277, NA, NA, NA, NA, 46.3063850402832, 45.5805702209473,
44.9196586608887, 42.2495613098145, 39.3051452636719, 38.7914810180664,
38.6069412231445, NA, NA, 44.6782455444336, 46.4024772644043,
44.4720573425293, 41.7361183166504, 42.3378067016602, 41.0018348693848,
39.3579216003418, 41.6303863525391, NA, 43.8207550048828,
46.0460357666016, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
3.32185006141663, NA, NA, NA, NA, NA, 4.98925733566284,
4.35699367523193,
4.94798421859741, 3.14325952529907, 2.93211793899536, NA,
NA, NA, NA, 4.52736520767212, 4.99723243713379, 5.13944292068481,
3.92965626716614, 3.43465113639832, NA, NA, NA, 3.55617475509644,
3.4659411907196, 5.24469566345215, 5.36995029449463, 4.61549234390259,
4.82002925872803, NA, NA, NA, NA, 4.20452928543091, 4.71502685546875,
5.20452785491943, 5.05676746368408, 5.9952244758606, 6.16778612136841,
4.69053316116333, NA, NA, 2.62325501441956, 4.74775457382202,
4.93133020401001, 5.02366256713867, 5.74016952514648, 6.28353786468506,
4.67424774169922, 4.56812858581543, NA, 1.88153350353241,
4.31531000137329, NA, NA, NA, NA, NA, NA), .Dim = c(63L,
4L)), offset = 0, gain = 1, inmemory = TRUE, fromdisk = FALSE,
nlayers = 4L, dropped = NULL, isfactor = c(FALSE, FALSE,
FALSE, FALSE), attributes = list(), haveminmax = TRUE,
min = c(12.7868213653564, 10.6471328735352, 38.6069412231445,
1.88153350353241), max = c(123.040382385254, 1148.87768554688,
46.4024772644043, 6.28353786468506), unit = "", names = c("ntl",
"pop", "tirs", "agbh")), legend = new(".RasterLegend",
type = character(0), values = logical(0), color = logical(0),
names = logical(0), colortable = logical(0)), title = character(0),
extent = new("Extent", xmin = 11878500, xmax = 11883000,
ymin = 1799000, ymax = 1802500), rotated = FALSE, rotation =
new(".Rotation",
geotrans = numeric(0), transfun = function ()
NULL), ncols = 9L, nrows = 7L, crs = new("CRS", projargs =
NA_character_),
srs = character(0), history = list(), z = list())
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Error when running multiscale GWR: Error in gw_weight_vec: Not compatible > with requested type: [type=NULL; target=double]
I am trying to run multiscale geographically weighted regression (MGWR) using the GWmodel package in R. When running the function gwr.multiscale this error is shown: Error in gw_weight_vec(vdist, bw, kernel, adaptive): Not compatible with requested type: [type=NULL; target=double]. An example: library(GWmodel) data(LondonHP) dist <- gw.dist(coordinates(londonhp)) ab_gwr <- gwr.multiscale(PURCHASE ~ FLOORSZ + PROF, data = londonhp, criterion = "dCVR", kernel = "gaussian", adaptive = FALSE, var.dMat.indx = 2, bws0 = c(100, 100, 100), bw.seled = rep(T, 3), dMats = list(dist, dist, dist), parallel.method = "omp", parallel.arg = "omp") I am following the example from the package's PDF. Also, I tried with my dataset and many different combinations (different number of covariates, the dMats parameter, adaptive bandwidths etc etc). -- Tziokas Nikolaos Cartographer Tel:(+44)07561120302 LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130> [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] GWmodel - create prediction raster at a finer spatial scale when using multiple independent variables
I am using the package GWmodel to predict at a finer spatial scale.
Basically, I want my output to be a single raster layer. For the GWR, I am
using one dependent and 2 independent variables. I don't know how to create
the prediction raster.
When using only one independent variable it's easy. For example:
library(GWmodel)
library(sp)
library(raster)
wd = "path/"
provoliko = "EPSG:7767"
# this is my fine res raster
tirs = raster(paste0(wd, "tirs.tif"))
regpoints <- as(tirs, "SpatialPoints")
# this df contains the coarse res variables
block.data = read.csv(paste0(wd, "block.data.csv"))
coordinates(block.data) <- c("x", "y")
proj4string(block.data) <- provoliko
eq1 <- ntl ~ tirs
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) <- provoliko
# slope
slope = as.data.frame(sf$tirs)
slope = SpatialPointsDataFrame(data = slope, coords = regpoints)
gridded(slope) <- TRUE
slope <- raster(slope)
raster::crs(slope) <- provoliko
gwr_pred = intercept + slope * s
writeRaster(gwr_pred,
paste0(wd, "ntl_gwr.tif"),
overwrite = TRUE)
When using multiple independent variables, I tried to stack them before I
convert them to spatialPoints, but then I don't know what to set in the
slope parameter.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] package spgwr: apply model parameters to a finer spatial scale
In order to apply *GWR*'s model parameters to a finer spatial scale using
the *spgwr *package:
1. calculate *GWR *at the coarse scale
2. apply step 1 again using the parameters* fit.points*, *predictions *and
*fittedGWRobject*.
The code:
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)
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)
Στις Σάβ 10 Δεκ 2022 στις 10:56 π.μ., ο/η Nikolaos Tziokas <
nikos.tzio...@gmail.com> έγραψε:
> I using the *R* package *spgwr *to perform geographically weighted
> regression (GWR). I want to apply the model parameters to a finer spatial
> scale but I am receiving this error: *Error in validObject(.Object):
> invalid class “SpatialPointsDataFrame” object: number of rows in data.frame
> and SpatialPoints don't match*.
>
> When I use another package for GWR, called *GWmodel*, I do not have this
> issue. For example using the *GWmodel*, I do:
>
> 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_
[R-sig-Geo] package spgwr: apply model parameters to a finer spatial scale
I using the *R* package *spgwr *to perform geographically weighted
regression (GWR). I want to apply the model parameters to a finer spatial
scale but I am receiving this error: *Error in validObject(.Object):
invalid class “SpatialPointsDataFrame” object: number of rows in data.frame
and SpatialPoints don't match*.
When I use another package for GWR, called *GWmodel*, I do not have this
issue. For example using the *GWmodel*, I do:
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)
How can I apply the GWR model parameters to a finer spatial scale, using
the spgwr package?
Here is the code, using the *spgwr *package:
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)
The fine resolution raster:
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=100 +y_0=100
+datum=WGS84 +units=m +no_defs')
The csv can be downloaded from here
<https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7?usp=sharing>
.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Apply geographically weighted regression's model parameters to a finer spatial scale
I have two raster layers, one coarse resolution and one fine resolution. My
goal is to extract GWR's coefficients (intercept and slope) and apply them
to my fine resolution raster.
I can do this easily when I perform simple linear regression. For example:
library(terra)
library(sp)
# focal terra
tirs = rast("path/tirs.tif") # fine res raster
ntl = rast("path/ntl.tif") # coarse res raster
# fill null values
tirs = focal(tirs,
w = 9,
fun = mean,
na.policy = "only",
na.rm = TRUE)
gf <- focalMat(tirs, 0.10*400, "Gauss", 11)
r_gf <- focal(tirs, w = gf, na.rm = TRUE)
r_gf = resample(r_gf, ntl, method = "bilinear")
s = c(ntl, r_gf)names(s) = c('ntl', 'r_gf')
model <- lm(formula = ntl ~ tirs, data = s)
# apply the lm coefficients to the fine res raster
lm_pred = model$coefficients[1] + model$coefficients[2] * tirs
But when I run GWR, the slope and intercept are not just two numbers (like
in linear model) but it's a range. For example, below are the results of
the GWR:
*Summary of GWR coefficient estimates*:
Min. 1st Qu. Median 3rd Qu. Max.
Intercept -1632.61196 -55.79680 -15.9968315.01596 1133.299
tirs20 -42.43020 0.43446 1.80026 3.75802 70.987
My question is how can extract GWR model parameters (intercept and slope)
and apply them to my fine resolution raster? In the end I would like to do
the same thing as I did with the linear model, that is, *GWR_intercept +
GWR_slope * fine resolution raster*.
Here is the code of GWR:
library(GWmodel)
library(raster)
block.data = read.csv(file = "path/block.data00.csv")
#create mararate df for the x & y coords
x = as.data.frame(block.data$x)
y = as.data.frame(block.data$y)
sint = as.matrix(cbind(x, y))
#convert the data to spatialPointsdf and then to spatialPixelsdf
coordinates(block.data) = c("x", "y")#gridded(block.data) <- TRUE
# specify a model equation
eq1 <- ntl ~ tirs
dist = GWmodel::gw.dist(dp.locat = sint, focus = 0, longlat = FALSE)
abw = bw.gwr(eq1,
data = block.data,
approach = "AIC",
kernel = "tricube",
adaptive = TRUE,
p = 2,
longlat = F,
dMat = dist,
parallel.method = "omp",
parallel.arg = "omp")
ab_gwr = gwr.basic(eq1,
data = block.data,
bw = abw,
kernel = "tricube",
adaptive = TRUE,
p = 2,
longlat = FALSE,
dMat = dist,
F123.test = FALSE,
cv = FALSE,
parallel.method = "omp",
parallel.arg = "omp")
ab_gwr
You can download the csv from here
<https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7?usp=sharing>.
The fine resolution raster I am using:
tirs = rast(ncols=407, nrows=342, nlyrs=1, xmin=509600, xmax=550300,
ymin=161800, ymax=196000, names=c('tirs'), crs='EPSG:27700')
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Downsample (aggregate) raster by a non-integer factor, using a Gaussian filter
I have a fine resolution raster (100m pixel size) and I want to downsample
it (aggregate) to 460m pixel size. The downsampling should be done using a
Gaussian filter (local function). The scale factor is 4.6, so non-integer.
The sigma (std) parameter of the Gaussian filter should be expressed in
cells (pixels) (also the Gaussian filter in this manner).
To put it differently I am trying to simulate coarse data as though they
were measured with a coarse point spread function (PSF). The PSF is assumed
to be a Gaussian filter.
So, in theory, for each coarse pixel I need to go to its center and
calculate the weights (from the PSF/Gaussian filter) needed for each fine
pixel surround it. To do this I need to apply a transfer function (TF;
e.g., Gaussian filter) to the fine data, but with a very large width.
The function aggregate from the terra package does not allow the fact
argument to be non-integer, so I thought to create a template raster (a
raster with no pixel values) from a coarse resolution raster at 460m
spatial scale that I have.
All in all, I need to create a custom function which downsamples a raster
using a Gaussian filter (units in pixels).
So far I created the template raster:
library(terra)
fr = rast("path/tirs.tif") # raster to be upscaled
cr = rast("path/ntl.tif")
# create an empty raster with the same dim, extent and crs as the coarse
res raster
(er <- rast(ext(cr), resolution=res(cr)))
crs(er) <- crs(cr)
The raster layers I am using:
fr = rast(ncols=216, nrows=417, nlyrs=1, xmin=582700, xmax=604300,
ymin=1005700, ymax=1047400, names=c('B10_median'), crs='EPSG:7767')
cr = rast(ncols=48, nrows=91, nlyrs=1, xmin=582360, xmax=604440,
ymin=1005560, ymax=1047420, names=c('avg_rad'), crs='EPSG:7767')
Any recommendations on how I can proceed?
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] How to check which model of an ensemble algorithm has been selected to perform regression?
I am using the R package machisplin
<https://github.com/jasonleebrown/machisplin> (it's not on CRAN) to
downscale a satellite image. According to the description of the package:
*The machisplin.mltps function simultaneously evaluates different
combinations of the six algorithms to predict the input data. During model
tuning, each algorithm is systematically weighted from 0-1 and the fit of
the ensembled model is evaluated. The best performing model is determined
through k-fold cross validation (k=10) and the model that has the lowest
residual sum of squares of test data is chosen. After determining the best
model algorithms and weights, a final model is created using the full
training dataset.*
My question is how can I check which model out of the 6 has been selected
for the downscaling? To put it differently, when I export the downscaled
image, I would like to know which algorithm (out of the 6) has been used to
perform the downscaling.
Here is the code:
library(MACHISPLIN)
library(raster)
library(gbm)
evi = raster("path/evi.tif") # covariate
ntl = raster("path/ntl_1600.tif") # raster to be downscaled
##convert one of the rasters to a point dataframe to sample. Use any
raster input.
ntl.points<-rasterToPoints(ntl,
fun = NULL,
spatial = FALSE)
##subset only the x and y data
ntl.points<- ntl.points[,1:2]
##Extract values to points from rasters
RAST_VAL<-data.frame(extract(ntl, ntl.points))
##merge sampled data to input
InInterp<-cbind(ntl.points, RAST_VAL)
#run an ensemble machine learning thin plate spline
interp.rast<-machisplin.mltps(int.values = InInterp,
covar.ras = evi,
smooth.outputs.only = T,
tps = T,
n.cores = 4)
#set negative values to 0
interp.rast[[1]]$final[interp.rast[[1]]$final <= 0] <- 0
writeRaster(interp.rast[[1]]$final,
filename = "path/ntl_splines.tif")
I vied all the output parameters (please refer to *Example 2* in the
package description) but I couldn't find anything relevant to my question.
I have posted a question
<https://github.com/jasonleebrown/machisplin/issues/10> on GitHub as well.
>From here
<https://drive.google.com/drive/folders/17SnegdgZ0192nqtP-95O55ZNPLngNMQH?usp=sharing>
you
can download my images.
--
Tziokas Nikolaos
Cartographer
Tel:(+44)07561120302
LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to upscale (reduce the spatial resolution) an image using a Gaussian filter in R?
I managed to solve the issue using the OpenImageR package and the function
down_sample_image.
library(raster)
library(OpenImageR)
r = raster("path/tirs.tif")
m = as.matrix(r)
psf = down_sample_image(m,
factor = 4.6, # zoom factor
gaussian_blur = T,
gauss_sigma = 0.5 * 100) # sigma * pixel size
e <- extent(r)
m2r <- raster(psf)
extent(m2r) <- e
raster::crs(m2r) <- "EPSG:7767"
res(m2r)
writeRaster(m2r, "path/tirs460.tif")
The function seems to do exactly what the creator of the area-to-point
regression Kriging, Qunming Wang, did with in his Matlab code (
https://github.com/qunmingwang/Code-for-S2-fusion/blob/master/PSF_template.m).
That is, he multiplied the width of the PSF with scale factor.
Στις Τετ, 7 Σεπ 2022, 15:46 ο χρήστης Greg Snow <538...@gmail.com> έγραψε:
> The `image_resize` function in the `magick` package may be a simple
> option. A gaussian filter is one of the options. This uses the
> ImageMagick program behind the scenes, but is a quick and easy
> interface within R for images (raster data).
>
> On Tue, Sep 6, 2022 at 10:46 AM Nikolaos Tziokas
> wrote:
> >
> > Hi,
> >
> > I want to resample a raster from 15m to 460m (pixel size) using a
> Gaussian
> > filter.
> >
> > *The goal*
> >
> > I am having a coarse image which I want to downscale. I also have a fine
> > resolution band to assist the downscaling. The downscaling method I am
> > using is called geographically weighted area-to-point regression Kriging
> > (GWATPRK). The method consists of two steps:
> >
> >1. GWR
> >2. area-to-point Kriging on the GWR's residuals
> >
> > In order to perform GWR using raster data, those needs to have the same
> > pixel size. This means that, my fine resolution image needs to be
> upscaled
> > to match the spatial resolution of the coarse band. This upscaling of the
> > fine band needs to be done using a Gaussian kernel with σ=0.5(i.e., the
> > PSF).
> >
> > How can I upscale (reduce the spatial resolution) a satellite image
> using a
> > Gaussian kernel (i.e., point spread function)?
> >
> > For reference, I am following the paper The effect of point spread
> function
> > on downscaling continua
> > <https://www.sciencedirect.com/science/article/pii/S092427162030229X>
> where
> > the authors at p.253 in Eq (9) mention:
> >
> > *the coarse image produced by upscaling the corresponding fine band k
> using
> > a PSF.*
> >
> > I googled how I can achieve that but unfortunately I couldn't find any
> > solution. So to do this, how can I use this Gaussian filter to change the
> > resolution (pixel size) of my image?
> >
> > Here is the image
> > <
> https://drive.google.com/drive/folders/18_1Kshb8WbT04gwOw4d_xhfQenULDXdB?usp=sharing
> >
> > I
> > am trying to convolve.
> >
> >
> > Many thanks
> >
> >
> > --
> > Tziokas Nikolaos
> > Cartographer
> >
> > Tel:(+44)07561120302
> > LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
> --
> Gregory (Greg) L. Snow Ph.D.
> 538...@gmail.com
>
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
