Re: [GRASS-user] v.kernel - Points weighted with attribute ?

[Micha Silver]

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: Hi Micha, thanks for your repply. For my own need, I think your solution is enough and I would like to try if you agree to send the R script. Le 22/10/19 à 12:53, Micha Silver a écrit : On 10/21/19 1:33 PM, Tais wrote: Hi all, I would like to compute a raster density map ("heat map") from vector points using v.kernel. The points represent house addresses. However, I would need to weight the points with a value stored in the attribute table (number of inhabitants). I saw on the module manual page that this functionality is a "known issue". I imagined a solution to by-passe the limitation: duplicate each point as many time as needed (according to the value stored in the attribute table) and use this new layer directly in v.kernel. However, I think it will definitely not be the most efficient way to do the trick. Do you have an alternative in mind ? Of course the best solution would be that someone in the development community would have the time and the kindness to implement this function directly in the module (I'm not skilled in C and thus cannot try myself). Maybe not the answer you are looking for, but you can create a weighted kernel density in R. I can send an R script that is intended to be run from within a GRASS session. It accepts a point vector and attrib column, and creates a kernel density raster, weighted by the attribute values (using the R defaults for bandwidth and Gaussian kernel) Requires, of course, that R is installed, with the packages: spatstat, raster, rgrass7. NB: Sorry if I should have posted this on the developer mailing list. I hesitated but decided to post here finally. Best, -- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918 -- 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