Wie wählt man die Rasterinterpolationsmethode?

Ich arbeite mit einem Datensatz von ~ 1,3 Millionen Wohngebäuden; Jede von ihnen wird durch eine normalverteilte Variable im Bereich von 0 bis 100 beschrieben. Der Datensatz deckt die gesamte Schweiz ab, sodass es offensichtlich Gebiete mit sehr dichter und sehr geringer Punktdichte gibt.

Die von unserem Datenanbieter auferlegte Einschränkung besteht darin, dass wir keine detaillierteren Datensätze als 5 Gebäude veröffentlichen können.

Ich habe darüber nachgedacht, diese Punktdaten in eine Rasterfläche umzuwandeln. Mein Hauptziel in dieser Übung wäre es, die Oberfläche zu erstellen, die den Fehler minimiert, der auftritt, wenn jemand versucht, den Punktwert aus diesem Raster zu ermitteln, indem er Daten an den interessierenden Punktpositionen seines Punktes extrahiert (im Vergleich zum Extrahieren aus dem Originaldatensatz von Punktpositionen).

Ich würde gerne wissen, welche Methode dafür am besten geeignet ist. In der Spatial Analyst- Toolbox stehen nur wenige zur Auswahl, in der geostatistischen Toolbox sogar noch mehr. Daher würde ich mich sehr über eine Hilfe bei den ersten Schritten freuen .

Es scheint, dass diese Frage mit einer früheren zusammenhängt , in der nach der Verschleierung solcher Daten mithilfe eines unregelmäßigen Rasters gefragt wird . Wenn wir akzeptieren, dass ein reguläres Raster verwendet wird, dann scheint es so

• Die meisten Zellen sollten groß genug sein, um fünf oder mehr Gebäude abzudecken

• Wenn Zellen nicht fünf Gebäude abdecken, sollten ihre Werte auf unvorhersehbare (aber kontrollierte) Weise geändert werden.

Wie Sie den Fehler messen, bestimmt die beste Lösung. Der in einer Zelle zu berechnende Wert sei y und die Werte der Gebäude in dieser Zelle (oder zumindest überlappend) seien x1 , x2 , ..., xk . Angenommen, jedes Gebäude hat ein nicht negatives "Interesse" (das proportional zur Anzahl der Bewohner sein kann), das als Ersatz für die erwartete Häufigkeit verwendet wird, mit der Ihr Raster geschätzt wird Gebäudewert. Nennen wir diese Ebenen w1 , ..., wk und w bezeichnen ihre Summe (ungleich Null).

1. Der durchschnittliche absolute Fehler ist das arithmetische Mittel der Fehlergrößen | y - xi | wie ich über die Gebäudeindizes reicht. Dies wird minimiert, indem y als Median des xi gewählt wird .

2. Der maximale Fehler ist der größte unter max (| y - xi |), da i über den Gebäudeindizes liegt. Dies wird minimiert, indem y als Mittelbereich von ( xi ) (Durchschnitt von max und min) gewählt wird. Dies wird jedoch stark von nur einem einzigen abweichenden Wert beeinflusst, sodass der Median möglicherweise vorzuziehen ist.

3. Der erwartete Fehler ist der gewichtete Durchschnitt von | y - xi | mit den Gewichten von wi / w . Dies wird minimiert, indem y als gewichteter Median des xi angenommen wird (aber kein GIS führt diese Berechnung für Sie durch - Sie müssen für solche Arbeiten ein statistisches oder mathematisches Paket wie Roder Mathematica verwenden .)

4. Der erwartete quadratische Fehler ist der gewichtete Durchschnitt von ( y - xi ) ^ 2. Es wird minimiert, indem y als gewichteter Mittelwert von xi angenommen wird , der gleich der Summe von wi xi / w ist .

Sie könnten mit (1) oder (2) aufgrund ihrer Einfachheit und direkten Interpretation zufrieden sein; Ich habe (3) und (4) eingefügt, um einen Eindruck von den Optionen zu bekommen. Um (1) zu implementieren, können Sie zunächst alle Daten mit einer Zellengröße gittern, die so klein ist, dass jedes Gebäude seine eigene Zelle belegt. (Bei einer Ausdehnung von etwa 200 mal 300 km würde eine Zellengröße von beispielsweise 5 m ein enormes Gitter von 40.000 mal 60.000 Zellen erfordern, aber nur etwa eine Million von ihnen wären belegt, wodurch nur etwa 10 MB Festplattenspeicher im nativen Zustand erforderlich wären Bogenformat, wenn Sie darauf achten, die Werte als Ganzzahlen zu speichern.) Aggregieren Sie dieses Raster mithilfe der Option zu einer größeren ZellengrößeMedianMöglichkeit. (Die Zellengröße des aggregierten Gitters würde wahrscheinlich etwa 100 m betragen, was ein landesweites Gitter von 2000 mal 3000 Zellen ergibt: ausreichend klein, um die nachstehend beschriebenen Verfahren nicht nur praktikabel, sondern auch schnell auszuführen.)

Sie sollten auch ein binäres Indikatorraster der Gebäude aggregieren - Sumdiesmal anfordern -, um die Anzahl der Gebäude pro Zelle zu zählen. Bei aggregierten Zellen mit einer Anzahl von weniger als 5 wird der Median zufällig gestört. Tun Sie dies mit einer ConOperation. Eine effektive, wenn auch etwas komplizierte Wahl für die Störung wäre, normalverteiltes Rauschen zum Logit des Werts hinzuzufügen (skaliert von 0 auf 1 anstatt von 0 auf 100): Dies garantiert ein Ergebnis, das immer noch zwischen 0 und 100 liegt Sie können auch alle Zellen leicht stören , sodass niemand gestörte Zellen von ungestörten Zellen unterscheiden kann, indem Sie die niedrigstwertigen Ziffern untersuchen.

Der Workflow für diese "logistische Störung" ist dann wie folgt. Es hängt von zwei Parametern ab: "Sigma" ist das Ausmaß der Störung der Zellen, die es benötigen, und "Epsilon" ist das minimale Ausmaß, um alle Zellen zu stören. Beides sind nicht negative Zahlen. Experimentieren Sie mit kleinen Teilgittern, beginnend mit Sigma = 0,15 und Epsilon = 0,01, und variieren Sie diese Parameter, bis die Ergebnisse zufriedenstellend sind. (Wenn Sie epsilon auf Null setzen, wird die Störung für solche Zellen insgesamt beseitigt.)

1. Beginnen Sie mit einem Raster [Z] von Medianwerten (alle im Bereich von 0 bis 100) und einem weiteren Raster [N], das die Anzahl der Gebäude in jeder Zelle zählt, die beide von erzeugt werden Aggregate.

2. Erstellen Sie mit einem ConBefehl wie ein Raster für die Störungsbeträge

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

3. Generieren Sie die normalverteilte Störung, indem Sie ein Raster mit Einheitsnormalvariablen verwenden (verwenden Sie CreateNormalRaster und multiplizieren Sie es mit dem vorherigen Raster. Nennen Sie das Ergebnis beispielsweise "e".

4. Berechnen Sie die gestörten Protokolle der Werte als

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

5. Zurück in Werte im Bereich 0..100 konvertieren:

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

Zur Veranschaulichung finden Sie hier einen RCode, mit dem Sie kleine Stichprobenaggregatgitter erstellen, stören und die gestörten mit den ursprünglichen Werten vergleichen können.

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)