Matching georeferenzierten Daten mit Shape-Datei in R

Question

Ich habe eine georeferenzierte Ereignisdatensatz in einem Datenrahmen der Form:

LONGITUDE LATITUDE VAR1 
33.4  4.4  5 
33.4  4.4  3 
33.4  4.4  1 
30.4  4.2  2 
28.4  5.1  2 

Es zählt Todesfälle in Ereignisse, die georeferenziert werden. Abgesehen davon, habe ich eine Formdatei der Provinzen in einem Land, wie folgt aus:

> str(shapefile) 
'data.frame': 216 obs. of 5 variables: 
$ CONSTI_COD: num 1 2 3 4 5 6 7 8 9 10 ... 
$ Area  : num 20 11.7 10.7 223.3 38.7 ... 
$ PROVINCE_NAME : Factor w/ 216 levels "CENTRAL","COAST",..: 4 4 4 4 4 4 4 4 2 2 ... 
$ Shape_Leng: num 0.193 0.152 0.201 0.872 0.441 ... 
$ Shape_Area: num 0.001628 0.000947 0.000867 0.018135 0.003145 ... 

..@ polygons :List of 216 
    .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots 
    .. .. .. ..@ Polygons :List of 1 
    .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots 
    .. .. .. .. .. .. ..@ labpt : num [1:2] 36.9 -1.3 
    .. .. .. .. .. .. ..@ area : num 0.00163 
    .. .. .. .. .. .. ..@ hole : logi FALSE 
    .. .. .. .. .. .. ..@ ringDir: int 1 
    .. .. .. .. .. .. ..@ coords : num [1:151, 1:2] 36.8 36.8 36.8 36.9 36.9 ... 
    .. .. .. ..@ plotOrder: int 1 
    .. .. .. ..@ labpt : num [1:2] 36.9 -1.3 
    .. .. .. ..@ ID  : chr "0" 
    .. .. .. ..@ area  : num 0.00163 
[...etc] 

Was ich tun müssen, ist die Ereignisdaten in den Provinzen zu platzieren, also fügen Sie eine vierte Spalte mit der ersten Daten Rahmen, der angibt, in welcher Provinz jedes Ereignis stattfand, basierend auf den Koordinaten. Also würde ich etwas in der Art haben:

LONGITUDE LATITUDE VAR1 PROVINCE 
33.4  4.4  5 CENTRAL 
33.4  4.4  3 CENTRAL 
33.4  4.4  1 CENTRAL 
30.4  4.2  2 COAST 
28.4  5.1  2 COAST 

Ist das möglich? Ich denke, ich habe vor einiger Zeit einen Beitrag gefunden, der erklärt, wie man das macht (außerhalb von Stack Overflow), aber ich kann es jetzt nicht finden.

Danke!

(Sorry, wenn es hier eine ähnliche Frage gibt. Ich habe eine Suche durchgeführt, aber ich habe keine Antwort gefunden, vielleicht weil ich nicht wirklich weiß, wonach ich suche. Ich würde einen Link sehr schätzen zu einem ähnlichen Beitrag.)

Answer 1

Worüber Sie sprechen ist eine "räumliche Verknüpfung" (oder "räumliche Schnittmenge" oder "Overlay"). Dies ist ziemlich einfach mit Hilfe der over Funktion aus dem sp Paket.

Hier ist ein Beispiel.

Zuerst laden und importieren wir ein Polygon Shapefile der Länder der Welt.

download.file(paste0('http://www.naturalearthdata.com/http//', 
        'www.naturalearthdata.com/download/110m/cultural/', 
        'ne_110m_admin_0_countries.zip'), 
       f <- tempfile()) 
unzip(f, exdir=tempdir()) 
library(rgdal) 
countries <- readOGR(tempdir(), 'ne_110m_admin_0_countries') 

Jetzt werden wir einige zufällige Koordinatendaten erstellen, die in das Ausmaß des Polygon Shapefile fallen. Wir definieren dann die Spalten x und y als coordinates und weisen dasselbe CRS wie das der Polygone zu (obwohl dies für Ihre Daten nicht der Fall sein kann; achten Sie darauf, korrekte Koordinatensysteme zuzuweisen).

pts <- data.frame(x=runif(10, -180, 180), y=runif(10, -90, 90), 
        VAR1=LETTERS[1:10]) 
coordinates(pts) <- ~x+y # pts needs to be a data.frame for this to work 
proj4string(pts) <- proj4string(countries) 

plot(countries) 
points(pts, pch=20, col='red') 

Jetzt können wir die räumliche Überlagerung durchführen:

over(pts, countries)$admin 

# [1] <NA>  <NA>  Turkey <NA>  
# [5] Macedonia <NA>  China  Argentina 
# [9] <NA>  Canada 
# 177 Levels: Afghanistan Albania ... Zimbabwe 

Beachten Sie, dass einige der zufälligen Punkten in diesem Fall in den Ozean gefallen (das heißt außerhalb Polygone). Wenn sie mit dem Polygonobjekt geschnitten werden, geben diese Punkte NA zurück.

Jetzt cbind wir das gewünschte Attribut zu pts:

cbind.data.frame(pts, country=over(pts, countries)$admin) 

#    x   y VAR1 country 
# 1 -52.59404 -37.422879 A  <NA> 
# 2 -33.88867 -40.194482 B  <NA> 
# 3 38.84383 37.272460 C Turkey 
# 4 -84.04949 7.118878 D  <NA> 
# 5 20.98272 40.920470 E Macedonia 
# 6 -155.32951 -37.612497 F  <NA> 
# 7 99.40166 38.630049 G  China 
# 8 -61.84025 -27.412885 H Argentina 
# 9 -37.65287 -3.666080 I  <NA> 
# 10 -112.81197 59.959475 J Canada