Generating raster (cost) surface with cell values being probabilities that generate weibull distribution

Question

I want to generate a raster surface (cost surface) based only on distance from center of extent where the values in the raster cells generate a Weibull distribution. I am interested specifically in the Weibull distribution, but the approach should be flexible enough to be modified to use with other distributions.

Here is some example code:

library(raster)

mat = matrix(data = rep(0,400),nrow=20,ncol=20)
rast = raster(mat,xmn=-10,xmx=10,ymn = -10,ymx=10)
crs(rast) = '+proj=utm +zone=12 +datum=WGS84'
plot(rast)

distFromPt = 1/raster::distanceFromPoints(rast, c(0,0)) 
plot(distrFromPt)

I want to be able to sample from the distrFromPt, such that the values of distrFromPt represent the probability that the cell would be chosen. For example, if the shape and scale parameters were both 1, then sampling 1000 random points should generate a distance histogram that looked like a weibull distribution.

hist(rweibull(1000,1,1),breaks=50)

I think the answer is somewhere in turning the values of distFromPt into quantiles, but just can't quite think how to do it.

Update:

Based on @Spacedman recommendations and another posting by @ Spacedman here, I have generated the following code:

library(raster)
mat = matrix(data = rep(0,400),nrow=20,ncol=20)
rast = raster(mat,xmn=-10,xmx=10,ymn = -10,ymx=10)
crs(rast) = '+proj=utm +zone=12 +datum=WGS84'
plot(rast)

distFromPt = 1/raster::distanceFromPoints(rast, c(0,0))
plot(distFromPt)

D = values(distFromPt)
Zw = qweibull(seq(0,1,len = length(rast)),1,1)
ZD = Zw[rank(D)]

values(rast) = ZD/sum(ZD)

hs = res(rast)/2

ptscell = sample(1:length(rast), 1000, prob=rast[], replace=TRUE)

centres = xyFromCell(rast,ptscell)

df = cbind(runif(nrow(centres),centres[,1]-hs[1],centres[,1]+hs[1]),
            runif(nrow(centres),centres[,2]-hs[2],centres[,2]+hs[2]))

pts = SpatialPointsDataFrame(coords = df,proj4string = crs(rast),
                               data = as.data.frame(df))

pts$Dist = pointDistance(pts,c(0,0),lonlat = FALSE, type = "Euclidean")

However, when you map the histogram of the distances, it doesn't look like I am able to recover a Weibull distribution.

hist(pts$Dist)

I think that this is because, even though the distant cells have a low probability of being selected, there are so many of them they "overwhelm" the cells close to the centroid that have a high probability, but low number of instances.

Answer 1

Suppose I have a vector D of 1000 elements from some distribution I don't know about:

> D = c(rnorm(500), rnorm(500,6,1))
> hist(D)

If I take 1000 points from quantiles of a Weibull(1,1) distribution:

> Zw = qweibull(seq(0,1,len=1000),1,1)
> hist(Zw)

Then I get the values you want. Note how the distribution of D has had no effect on the output distribution.

I suspect the next step you want is to match the values in your raster to values from that Weibull by ranking them so they are ordered, so you have a monotonic transformation from D to Zw.

> ZD = Zw[rank(D)]
> plot(D, ZD)

and this is just an empirical QQ plot of D and ZD. If you knew the distribution of D then you could work out the theoretical QQ plot and use that as a function to compute ZD for each D value, but maybe you don't need that.

What we've sort of done here is first totally flatten the distribution of the data by ranking it 1:N, then taking N values from the target distribution by quantiles (not by random sampling). So this can be done for any source or target distribution (with some diddling for discrete distributions...).