Untitled

Benötigten libraries und daten importieren:

library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(ggplot2) # optional, nur für Visualisierung
punkte_a <- read_sf("data/Punkte_A.shp")
punkte_b <- read_sf("data/Punkte_B.shp")

border = st_read("data/Grenze.shp")
## Reading layer `Grenze' from data source `/home/rstudio/data/Grenze.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -558211.2 ymin: 282997.4 xmax: -557225.3 ymax: 283903.5
## Projected CRS: NAD_1983_Albers

Daten vorbereiten:

st_crs(punkte_a) == st_crs(punkte_b)
## [1] FALSE
st_crs(punkte_a) == st_crs(border)
## [1] TRUE
st_crs(punkte_b) == st_crs(border)
## [1] FALSE

# alle datensätze müssen im gleichen CRS sein
# ich konvertiere punkte_b ins gleiche crs wie punkte_a
punkte_b <- st_transform(punkte_b, st_crs(punkte_a))

st_crs(punkte_a) == st_crs(punkte_b)
## [1] TRUE
st_crs(punkte_b) == st_crs(border)
## [1] TRUE

Kurz schauen wie’s aussieht:

p <- ggplot(data = punkte_a) + 
  geom_sf(data = border) +
  geom_sf(data = punkte_b, aes(col = "punkte_b"), pch = 3) +
  geom_sf(aes(col = "punkte_a"))  +
  geom_sf_text(aes(label = Plot_Numbe), nudge_x = 15, nudge_y = 15)

p

Frage 1

Ich möchte die Distanzen zwischen den Punkten in A und B berechnen. Das Ziel am Schluss ist den nächsten Punkt aus A zu jedem einzelnen Punkt aus B zu finden.

Wenn ich dich richtig verstehe, willst du für jeden Punkt in B den nächsten aus A finden? Wenn das stimmt ist die Antwort ziemlich langweilig: der punkt mit der Plot_Numbe 4 ist der nächste Punkt zu allen Punkten B. Das ergibt sich aus der obigen Visualisierung (dieser Punkt liegt in mitten aller punkte_b.

Berechnen kannst du das folgendermassen:

st_nearest_feature(punkte_b, punkte_a)
##  [1] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [39] 4

Wenn ich dich falsch verstanden habe und du das ganze umgekehrt meinst, musst du nur die Reihenfolge umdrehen.

st_nearest_feature(punkte_a, punkte_b)
##  [1] 24 24 25 17 39 14 14 14 14 14 14 14 14 14 20 23 14 14 20 23 23 36 36 36 36
## [26] 23 23 23 23 39 20 14 24 36 14 14  3 36

Du kannst das ganze auch etwas “händischer” berechnen und hättest so mehr kontrolle über die Resultate. (Zum Beispiel könntest du auch die zweitnächsten Features ermitteln oder weisst auch die jeweiligen Distanzen)

dist_matrix <- st_distance(punkte_a, punkte_b)

# Die resultierende Matrix hat soviele `rows` wie `punkte_a` features hat und so viele `columns` wie `punkte_b` features hat. 
nrow(punkte_a)
## [1] 38
nrow(dist_matrix)
## [1] 38

nrow(punkte_b)
## [1] 39
ncol(dist_matrix)
## [1] 39

# Nun kannst du von der obigen matrix *pro Spalte* schauen, welche Zeile den kleinsten Wert aufweist. 
apply(dist_matrix, MARGIN = 2, FUN = which.min)
##  [1] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [39] 4

# ist übrigens das gleiche wie:
st_nearest_feature(punkte_b, punkte_a)
##  [1] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
## [39] 4
# Margin = 1 geht durch die Zeilen von dist_matrix. Margin = 2 geht über die Spalten
# which.min gibt die *Position* des kleinsten wertes zurück. "min" würde dir lediglich den Wert zurück geben.

# Oder du schaust es pro Zeile an
apply(dist_matrix, MARGIN = 1, FUN = which.min)
##  [1] 24 24 25 17 39 14 14 14 14 14 14 14 14 14 20 23 14 14 20 23 23 36 36 36 36
## [26] 23 23 23 23 39 20 14 24 36 14 14  3 36

# ist übrigens das gleiche wie:
st_nearest_feature(punkte_a, punkte_b)
##  [1] 24 24 25 17 39 14 14 14 14 14 14 14 14 14 20 23 14 14 20 23 23 36 36 36 36
## [26] 23 23 23 23 39 20 14 24 36 14 14  3 36

Das ganze etwas grafischer:

# Für jeden Punkt in B, welches ist der nächste aus A?
p +
  geom_sf(data = st_nearest_points(punkte_b, st_union(punkte_a)))

# Für jeden Punkt in A, welches ist der nächste aus B?
p +
  geom_sf(data = st_nearest_points(st_union(punkte_b), punkte_a))

Frage 2

Ich habe basiert auf einem Shapefile A thiessenpolygone/voronoipolygone generiert. Dafür habe ich den Code aus dem ReMe Tutorial genutzt:

punkte_a_union <- st_union(punkte_a)

thiessenpolygone <- st_voronoi(punkte_a_union)

thiessenpolygone <- st_cast(thiessenpolygone)

thiessenpolygone_clip <- st_intersection(thiessenpolygone,border)

thiessenpolygone_clip <- st_sf(thiessenpolygone_clip)
thiessenpolygone_clip$id <- seq_len(nrow(thiessenpolygone_clip))

ggplot(thiessenpolygone_clip) +
  geom_sf() +
  geom_sf(data = punkte_b, alpha = .2) +
  geom_sf_text(aes(label = id), colour = "red") 

Nun möchte ich prüfen ob punkte aus B auf einem dieser Polygone liegen. Das Ziel wäre, dass ich nachher für jeden Eintrag in B die Nummer des jeweiligen Polygon steht. Dafür wollte ich folgenden Code nutzen:

trees_UAV_with_plot_nr = st_intersection(thiessenpolygone_clip , punkte_b) 
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

st_intersection macht eine intersection der beiden Geometrien. Was du aber brauchst ist eine Antwort auf die Frage: welches Polygon intersected mit dem jeweiligen Punkt. Das ist nicht st_intersection sondern st_intersects.

punkte_b_intersects <- st_intersects(punkte_b, thiessenpolygone_clip)

st_intersects retourniert eine Sparse geometry, also eine Liste mit der gleichen Länge wie punkte_b. Jeder eintrag in der Liste ist ein nummerischer Wert und zeigt an, welches Feature aus thiessenpolygone_clip mit dem entsprechend punkt_b interesected.

In deinem Fall is dies immer 1 Feature, nicht mehr und nicht weniger (in anderen Situationen könnten es auch 0 oder ganz viele sein). Aber in deinem Fall kannst du die Liste mit unlist() in einen Zahlenvektor konvertieren und auch gerade als neue Spalte zu punket_b hinzufügen:

unlist(punkte_b_intersects)
##  [1] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
## [39] 7

punkte_b$thiessen_pos <- unlist(punkte_b_intersects)

Jetzt sind die Polygone jedoch eine Liste und nicht einfach ein “Geometriedatensatz”, und lässt sich nicht mit den Punkten verschneiden.

Geometrien sind in R immer als Liste abgespeichert.

is.list(punkte_a$geometry)
## [1] TRUE
is.list(punkte_b$geometry)
## [1] TRUE
is.list(border$geometry)
## [1] TRUE