I am trying to optimize a conductance surface using multiple categorical and one (hsm_ag) continuous layers. I noticed that the 'surface' object that is produced by the conductance_surface() function has all NA's except for the "elevation" column. If I remove the "elevation" layer from the RasterStack before feeding it to the conductance_surface() function, the resultant 'surface' object contains categorical values, as expected. Then the radish function seems to be satisfied with this. Is there any idea why I am getting NAs for everything with certain landscape layers included? Below is my code.
ev_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/elev_raw_agg5.tif")
lu_ag <-raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/landuse_raw_agg5.tif")
rd_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/roads_raw_agg5.tif")
rio_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/rios_raw_agg5.tif")
cit_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/cities_raw_agg5.tif")
camp_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/camp_raw_agg5.tif")
hsm_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/forest_raw_agg5.tif")
fc_ag <-raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/HSM_RAW_FINAL_agg5.tif")
make_factor <- function(r){
if (!is.factor(r)) {
r <- calc(r, function(x) { x[!is.na(x)] <- as.integer(round(x[!is.na(x)])); x })
r <- ratify(r)
}
Optional: name the levels via VALUE column (if you have a mapping)
RAT <- levels(r)[[1]]
RAT$VALUE <- paste0("class_", RAT$ID) # <-- replace with your real labels
levels(r) <- RAT
r
}
lu_ag <- make_factor(lu_ag)
rd_ag <- make_factor(rd_ag)
rio_ag <- make_factor(rio_ag)
cit_ag <- make_factor(cit_ag)
camp_ag<- make_factor(camp_ag)
hsm_ag <- make_factor(hsm_ag)
fc_ag <- make_factor(fc_ag)
Scale the continuous layer
ev_sc <- raster::scale(ev_ag)
Build a RasterStack mixing factors + continuous
covariates <- raster::stack(list(
elevation = ev_sc,
landuse = lu_ag,
roads = rd_ag,
rivers = rio_ag,
cities = cit_ag,
campesinas = camp_ag,
hsm = hsm_ag,
forestcover = fc_ag
))
Force an identical NA mask across ALL layers (radish warns if not)
mask_common <- raster::calc(covariates, function(v) if (all(!is.na(v))) 1L else NA_integer_)
covariates <- raster::mask(covariates, mask_common)
covariates <- raster::stack(covariates) # ensure class is exactly RasterStack
make sites spatial points
sites_pts <- as(sites_sp, "SpatialPoints")
if (!compareCRS(covariates, sites_pts)) {
sites_pts <- spTransform(sites_pts, crs(covariates))
}
cells <- cellFromXY(covariates[[1]], coordinates(sites_pts))
keep <- !is.na(cells) & !is.na(raster::extract(covariates[[1]], sites_pts))
if (any(!keep)) message(sum(!keep), " site(s) fall on NA/outside; dropping.")
sites_pts <- sites_pts[keep, ]
Build the conductance surface
surface <- conductance_surface(
covariates = covariates, # RasterStack with factors + scaled continuous
coords = sites_pts, # SpatialPoints
directions = 8
)
save output
outdir <- "/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/EditedLayers2Sept25/radish_out"
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
saveRDS(surface, file.path(outdir, "surface.rds"))
I am trying to optimize a conductance surface using multiple categorical and one (hsm_ag) continuous layers. I noticed that the 'surface' object that is produced by the conductance_surface() function has all NA's except for the "elevation" column. If I remove the "elevation" layer from the RasterStack before feeding it to the conductance_surface() function, the resultant 'surface' object contains categorical values, as expected. Then the radish function seems to be satisfied with this. Is there any idea why I am getting NAs for everything with certain landscape layers included? Below is my code.
ev_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/elev_raw_agg5.tif")
lu_ag <-raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/landuse_raw_agg5.tif")
rd_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/roads_raw_agg5.tif")
rio_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/rios_raw_agg5.tif")
cit_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/cities_raw_agg5.tif")
camp_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/camp_raw_agg5.tif")
hsm_ag <- raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/forest_raw_agg5.tif")
fc_ag <-raster("/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/landscape_data_north/raw/HSM_RAW_FINAL_agg5.tif")
make_factor <- function(r){
if (!is.factor(r)) {
r <- calc(r, function(x) { x[!is.na(x)] <- as.integer(round(x[!is.na(x)])); x })
r <- ratify(r)
}
Optional: name the levels via VALUE column (if you have a mapping)
RAT <- levels(r)[[1]]
RAT$VALUE <- paste0("class_", RAT$ID) # <-- replace with your real labels
levels(r) <- RAT
r
}
lu_ag <- make_factor(lu_ag)
rd_ag <- make_factor(rd_ag)
rio_ag <- make_factor(rio_ag)
cit_ag <- make_factor(cit_ag)
camp_ag<- make_factor(camp_ag)
hsm_ag <- make_factor(hsm_ag)
fc_ag <- make_factor(fc_ag)
Scale the continuous layer
ev_sc <- raster::scale(ev_ag)
Build a RasterStack mixing factors + continuous
covariates <- raster::stack(list(
elevation = ev_sc,
landuse = lu_ag,
roads = rd_ag,
rivers = rio_ag,
cities = cit_ag,
campesinas = camp_ag,
hsm = hsm_ag,
forestcover = fc_ag
))
Force an identical NA mask across ALL layers (radish warns if not)
mask_common <- raster::calc(covariates, function(v) if (all(!is.na(v))) 1L else NA_integer_)
covariates <- raster::mask(covariates, mask_common)
covariates <- raster::stack(covariates) # ensure class is exactly RasterStack
make sites spatial points
sites_pts <- as(sites_sp, "SpatialPoints")
if (!compareCRS(covariates, sites_pts)) {
sites_pts <- spTransform(sites_pts, crs(covariates))
}
cells <- cellFromXY(covariates[[1]], coordinates(sites_pts))
keep <- !is.na(cells) & !is.na(raster::extract(covariates[[1]], sites_pts))
if (any(!keep)) message(sum(!keep), " site(s) fall on NA/outside; dropping.")
sites_pts <- sites_pts[keep, ]
Build the conductance surface
surface <- conductance_surface(
covariates = covariates, # RasterStack with factors + scaled continuous
coords = sites_pts, # SpatialPoints
directions = 8
)
save output
outdir <- "/projectnb/vervpop/mazarate/MEL_LANDSCAPE_GEN/EditedLayers2Sept25/radish_out"
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
saveRDS(surface, file.path(outdir, "surface.rds"))