2013-03-13 20 views
5

Qualcuno sa come creare un dataset di elevazione raster di simulazione, ovvero una matrice 2d con valori di elevazione realistici, in R? R jitter non sembra appropriato. In Java/Elaborazione funzione noise() realizza questo con un algoritmo Perlin noise es .:Dati altimetrici simulati realistici in rumore R/Perlin

size(200, 200); 
float ns = 0.03; // for scaling 
for (float i=0; i<200; i++) { 
    for (float j=0; j<200; j++) { 
    stroke(noise(i*ns, j*ns) * 255); 
    point(i, j); 
    } 
} 

enter image description here

Ma ho trovato alcun riferimento al rumore Perlin in letteratura R. Grazie in anticipo.

+1

vedere né il pacchetto 'RandomFields' (metodi per la simulazione di una vasta varietà di corse gaussiane campi dom), o forse una superficie frattale: http://www.oceanographerschoice.com/2010/10/fractal-landscapes-in-r-part-two/ (Ho un vecchio codice per questo) –

+0

Una ricerca rapida [ richiama qualche codice C++] (http://www.dreamincode.net/forums/topic/66480-perlin-noise/), che non dovrebbe essere troppo difficile da regolare. Potresti usare 'Rcpp' per usarlo direttamente o tradurlo in R. – Roland

+0

Grazie Ben, Roland. Il blog oceanograderschoice è una cosa utile da provare ad un certo punto quando ho tempo, come se stessimo imparando C++. – geotheory

risposta

10

Ecco un'implementazione in R, seguendo le spiegazioni in http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html

perlin_noise <- function( 
    n = 5, m = 7, # Size of the grid for the vector field 
    N = 100, M = 100 # Dimension of the image 
) { 
    # For each point on this n*m grid, choose a unit 1 vector 
    vector_field <- apply(
    array(rnorm(2 * n * m), dim = c(2,n,m)), 
    2:3, 
    function(u) u/sqrt(sum(u^2)) 
) 
    f <- function(x,y) { 
    # Find the grid cell in which the point (x,y) is 
    i <- floor(x) 
    j <- floor(y) 
    stopifnot(i >= 1 || j >= 1 || i < n || j < m) 
    # The 4 vectors, from the vector field, at the vertices of the square 
    v1 <- vector_field[,i,j] 
    v2 <- vector_field[,i+1,j] 
    v3 <- vector_field[,i,j+1] 
    v4 <- vector_field[,i+1,j+1] 
    # Vectors from the point to the vertices 
    u1 <- c(x,y) - c(i,j) 
    u2 <- c(x,y) - c(i+1,j) 
    u3 <- c(x,y) - c(i,j+1) 
    u4 <- c(x,y) - c(i+1,j+1) 
    # Scalar products 
    a1 <- sum(v1 * u1) 
    a2 <- sum(v2 * u2) 
    a3 <- sum(v3 * u3) 
    a4 <- sum(v4 * u4) 
    # Weighted average of the scalar products 
    s <- function(p) 3 * p^2 - 2 * p^3 
    p <- s(x - i) 
    q <- s(y - j) 
    b1 <- (1-p)*a1 + p*a2 
    b2 <- (1-p)*a3 + p*a4 
    (1-q) * b1 + q * b2 
    } 
    xs <- seq(from = 1, to = n, length = N+1)[-(N+1)] 
    ys <- seq(from = 1, to = m, length = M+1)[-(M+1)] 
    outer(xs, ys, Vectorize(f)) 
} 

image(perlin_noise()) 

Perlin noise

Si può avere una struttura più frattale aggiungendo tali matrici, con differenti dimensioni della griglia.

a <- .6 
k <- 8 
m <- perlin_noise(2,2,2^k,2^k) 
for(i in 2:k) 
    m <- m + a^i * perlin_noise(2^i,2^i,2^k,2^k) 
image(m) 
m[] <- rank(m) # Histogram equalization 
image(m) 

Perlin noise, Sum with different grid sizes

+0

Spot su Vincent :) Ho passato un paio d'ore a guado attraverso altre implementazioni - questo codice è molto più elegante di quanto sarebbe successo al mio! Sai chi merita l'attribuzione se necessario? – geotheory

+0

Grazie anche per la modifica. Sono incuriosito dai ripetuti salti in cui _high_ e _low_ punti estremi si avvicinano l'un l'altro. Il modello sembra perfetto oltre a questi punti. – geotheory

+2

Questo artefatto era dovuto a un 'sqrt 'mancante durante la normalizzazione dei vettori; L'ho riparato e aggiornato di conseguenza le trame. –

0

Un metodo alternativo:

require(geoR) 
sim <- grf(441, grid="reg", cov.pars=c(1, .25)) 
image(sim, col=gray(seq(1, .1, l=30))) 

enter image description here

può estrarre i dati degli oggetti con cbind(sim[[1]], z = sim[[2]])