Como escolher o método de interpolação raster?

Eu trabalho com um conjunto de dados de ~ 1,3 milhão de edifícios residenciais; cada um deles é descrito por uma variável normalmente distribuída, variando de 0 a 100. O conjunto de dados abrange toda a Suíça, de modo que obviamente existem áreas com densidade de pontos muito densa e muito esparsa.

A limitação imposta pelo nosso provedor de dados é que não podemos liberar conjuntos de dados mais detalhados que cinco edifícios.

Eu estava pensando em converter esses dados de ponto em uma superfície raster. Meu principal objetivo neste exercício seria criar a superfície que minimize o erro que surgirá quando alguém tentar obter o valor em pontos dessa varredura extraindo dados nos locais de interesse de seus pontos (em comparação com a extração do conjunto de dados original de localizações dos pontos).

Gostaria de saber que tipo de método deve ser mais adequado para fazer isso. Existem poucas opções para escolher na caixa de ferramentas do Spatial Analyst e ainda mais na geoestatística, por isso gostaria muito de receber ajuda para começar.

Parece que esta pergunta está relacionada a uma pergunta anterior que pergunta sobre disfarçar esses dados usando uma grade irregular. Se aceitarmos que uma grade regular será usada, parece que

• A maioria das células deve ser grande o suficiente para cobrir cinco ou mais prédios e

• Quando as células não cobrem cinco edifícios, seus valores devem ser alterados de maneiras imprevisíveis (mas controladas).

Como você mede o erro determina a melhor solução. Seja o valor a ser computado em uma célula seja y e os valores dos prédios localizados dentro (ou pelo menos sobrepostos) dessa célula sejam x1 , x2 , ..., xk . Além disso, suponha que cada edifício tenha um "nível de interesse" não negativo (que pode ser proporcional ao número de moradores), para ser usado como substituto do número esperado de vezes que sua varredura será usada para estimar que valor do edifício. Vamos chamar estes níveis w1 , ..., wk , respectivamente, e deixe w designar sua soma (diferente de zero).

1. O erro absoluto médio é a média aritmética dos tamanhos dos erros | y - xi | como i varia sobre os índices de construção. Isso é minimizado escolhendo y como a mediana do xi .

2. O erro máximo é o maior entre max (| y - xi |), conforme i varia sobre os índices de construção. Isso é minimizado escolhendo y como a faixa intermediária de ( xi ) (média do máximo e do mínimo). Isso é fortemente influenciado por apenas um único valor externo, portanto, a mediana pode ser preferível.

3. O erro esperado é a média ponderada de | y - xi | com os pesos dados por wi / w . Esta é minimizado tomando y a ser uma média ponderada do xi (mas sem GIS vai fazer este cálculo para você - você terá que usar um pacote estatístico ou matemática como Rou Mathematica . Para esse trabalho)

4. O erro quadrático esperado é a média ponderada de ( y - xi ) ^ 2. É minimizado considerando y como a média ponderada do xi , igual à soma de wi xi / w .

Você pode ser feliz com (1) ou (2) devido à sua simplicidade e interpretação direta; Incluí (3) e (4) para dar uma idéia das opções. Para implementar (1), você pode começar agrupando todos os dados com um tamanho de célula tão pequeno que cada edifício ocupa sua própria célula. (Em uma extensão de cerca de 200 por 300 km, um tamanho de célula de, digamos, 5 m exigiria uma enorme grade de 40.000 a 60.000 células, mas apenas cerca de um milhão delas seria ocupada, exigindo, portanto, apenas cerca de 10 MB de armazenamento em disco no ambiente nativo. formato de arco se você tomar cuidado para armazenar os valores como inteiros.) agregar esta grade a um cellsize maior usando oMedianopção. (O tamanho da célula da grade agregada provavelmente seria de cerca de 100m, fornecendo uma grade de 2000 por 3000 células em todo o país: suficientemente pequena para tornar os procedimentos descritos abaixo não apenas praticáveis, mas de execução rápida.)

Você também deve agregar uma grade de indicadores binários dos edifícios - solicitando Sumesse tempo - para contar o número de edifícios por célula. Para células agregadas com contagens inferiores a 5, perturbe aleatoriamente a mediana. Faça isso com uma Conoperação. Uma escolha eficaz, embora um pouco complicada, para a perturbação seria adicionar ruído distribuído normalmente ao logit do valor (dimensionado de 0 a 1 em vez de 0 a 100): isso garantirá um resultado que ainda esteja entre 0 e 100 Você também pode perturbar levemente todas as células, para que ninguém possa distinguir células perturbadas de células imperturbáveis, inspecionando os dígitos menos significativos.

O fluxo de trabalho para esse procedimento de "perturbação logística" é o seguinte. Depende de dois parâmetros: "sigma" é a quantidade de perturbação das células que precisa e "epsilon" é a quantidade mínima para perturbar todas as células. Ambos são números não negativos. Experimente em pequenas subgrelhas, começando com sigma = 0,15 e epsilon = 0,01, e varie esses parâmetros até que os resultados sejam satisfatórios. (Definir epsilon para zero eliminará completamente a perturbação dessas células.)

1. Comece com uma grade [Z] de valores medianos (todos no intervalo de 0 a 100) e outra grade [N] contando o número de edifícios em cada célula, ambos produzidos por Aggregate.

2. Crie uma grade para as quantidades de perturbação usando um Concomando como

   Con["N" < 5, sigma, epsilon]
   

3. Gere a perturbação normalmente distribuída tomando uma grade de variáveis ​​normais da unidade (use CreateNormalRaster e multiplique pela grade anterior. Chame o resultado "e", por exemplo.

4. Calcule os logits perturbados dos valores como

   [Logit] = log("z" / (100 - "z")) + "e"
   

5. Converta de volta para valores no intervalo 0..100:

   100 / (1 + exp(-"logit"))
   

Para ilustrar, aqui está um Rcódigo para criar pequenas grades agregadas de amostra, perturbá-las e comparar os perturbados com os valores originais.

ncol <- 30; nrow <- 20
seed.random <- 17
x <- rpois(ncol * nrow, 5)
y <- floor(100 / (1 + exp(-(rnorm(ncol * nrow, mean = -2, sd = 1/sqrt(x))))))

sigma <- 0.15
epsilon <- 0.01

e <- rnorm(ncol*nrow, sd = ((x < 5)*sigma + (x >= 5)*epsilon))
logit <- log(y / (100 - y)) + e
y0 <- 100 / (1 + exp(-logit))

library(raster)
z <- matrix(y, ncol=ncol)
n <- matrix(x, ncol=ncol)
z0 <- matrix(y0, ncol=ncol)

par(mfrow=c(2,2))
n.r <- raster(n)
plot(n.r, main="Counts of residences [N]")

z.r <- raster(z)
plot(z.r, main="Median values [Z]")

z0.r <- raster(z0)
plot(z0.r, main="Perturbed median values")

plot(y, y0, type="n", xlab="Original medians", ylab="Perturbed medians",
     main="Perturbed vs. original medians")
points(y[x < 5], y0[x < 5], col="Red")
points(y[x >= 5], y0[x >= 5], pch=19)