dati georeferenziati abbinamento con file di forma in R

Question

Ho una serie di dati georeferenziata evento, in una cornice di dati del modulo:

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 

Conta decessi in eventi che sono georeferenziati. A prescindere da esso, ho un file di forma delle province in un paese, come questo:

> 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] 

Che cosa devo fare è di inserire i dati degli eventi all'interno delle province, cioè aggiungere una quarta colonna per i primi dati frame che indica in quale provincia è accaduto ogni evento, in base alle coordinate. Quindi avrei qualcosa di simile:

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 

È possibile? Penso di aver trovato qualche tempo fa un post che spiegava come fare questo (fuori dallo Stack Overflow), ma non riesco a trovarlo ora.

Grazie!

(Scusa se c'è una domanda simile qui. Ho fatto una ricerca, ma non ho trovato una risposta, forse perché non so davvero cosa sto cercando. Gradirei davvero un link a un post simile.)

Answer 1

Quello di cui stai parlando è un "join spaziale" (o "intersezione spaziale" o "sovrapposizione"). Questo è abbastanza semplice con l'aiuto della funzione over dal pacchetto sp.

Ecco un esempio.

Per prima cosa, scarichiamo e importiamo uno shapefile poligono dei paesi del mondo.

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') 

Ora creeremo alcuni dati di coordinate casuali che rientrano nell'estensione dello shapefile del poligono. Quindi definiamo le colonne x e come coordinates e assegniamo lo stesso CRS di quello dei poligoni (anche se questo potrebbe non essere il caso per i tuoi dati, assicurati di assegnare i sistemi di coordinate corretti).

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') 

Ora siamo in grado di eseguire la sovrapposizione spaziale:

over(pts, countries)$admin 

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

nota che in questo caso, alcuni dei punti casuali caddero nell'oceano (cioè poligoni esterni). Quando si intersecano con l'oggetto poligono, questi punti restituiscono NA.

Ora abbiamo cbind l'attributo desiderato per 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