Here is the script. Let me know how it goes. (I might try to take time to make it into a GRASS addon)
On 10/22/19 2:33 PM, Taïs wrote:
-- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918 |
# Weighted Kernel Density in R # execute from a GRASS session # Command line arguments: GRASS point vector, column of weights ##################################################################### suppressPackageStartupMessages(c(library(rgrass7), library(spatstat), library(sp), library(raster))) # Test case, in GRASS location nc_basic_spm, using the schools vector # This vector has attribut CORECAPACI, used as weight in this case args = commandArgs(trailingOnly=TRUE) if (length(args) < 2) { print("Enter GRASS point vector and weight attribute column as command line arguments") print("Exiting...") quit(save="no") } else { input_vect = args[1] output_rast = paste0(input_vect, "_density") weights_col = args[2] } use_sp() points = readVECT(input_vect) # Check geometry type and weight attribute column if (class(points)!= "SpatialPointsDataFrame") { print("Wrong geometry type. Must be POINT") quit(save="no") } weights_idx = which(names(points) == weights_col) if (is.empty(weights_idx)) { print(paste("No attribute column:", weights_col)) quit(save="no") } # Remove NA points = points[!is.na(points[[weights_idx]]),] wts = points[[weights_idx]] print(paste("Using:", length(wts), "points")) # Prepare ppp object pts_coords = coordinates(points) pts_ext = c(min(pts_coords[,1]), max(pts_coords[,1]), min(pts_coords[,2]), max(pts_coords[,2])) pts_ppp = as.ppp(pts_coords, pts_ext) # Density (default Gaussian kernel) pts_dens = density(pts_ppp, weights = wts) # Convert to SGDF to save back to GRASS pts_sgdf = as(raster(pts_dens), "SpatialGridDataFrame") writeRAST(x = pts_sgdf, vname = output_rast)
_______________________________________________ grass-user mailing list grass-user@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-user