Come eseguire l'analisi spaziale in R con sf

Dove voti? Chi siete legislatori? Qual è il tuo codice postale? Queste domande hanno qualcosa in comune dal punto di vista geospaziale: la risposta implica determinare in quale poligono cade un punto.

Tali calcoli vengono spesso eseguiti con software GIS specializzato. Ma sono anche facili da fare in R. Hai bisogno di tre cose:

  1. Un modo per geocodificare gli indirizzi per trovare latitudine e longitudine; 
  2. Shapefile che delineano i confini del poligono del codice postale; e 
  3. Il pacchetto sf.

Per la geocodifica, di solito utilizzo l'API geocod.io. È gratuito per 2.500 ricerche al giorno e ha un bel pacchetto R, ma per usarlo è necessaria una chiave API (gratuita). Per aggirare quel po 'di complessità per questo articolo, userò l'API Nominatim Open Street Map gratuita e open source. Non richiede una chiave. Il pacchetto tmaptools ha una funzione,, geocode_OSM()per usare quell'API .

Importazione e preparazione dei dati geospaziali

Userò i pacchetti sf, tmaptools, tmap e dplyr. Se vuoi continuare, carica ciascuno con pacman::p_load()o installa quelli non ancora sul tuo sistema con install.packages(), quindi carica ciascuno con library().

Per questo esempio, creerò un vettore con due indirizzi, il nostro ufficio a Framingham, Massachusetts, e l'ufficio RStudio a Boston.

indirizzi <- c ("492 Old Connecticut Path, Framingham, MA",

"250 Northern Ave., Boston, MA")

La geocodifica è semplice con geocode_OSM. Puoi vedere i risultati stampando le prime tre colonne incluse latitudine e longitudine:

geocoded_addresses <- geocode_OSM (indirizzi)

print (geocoded_addresses [, 1: 3])

query lat lon

# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105

# 2250 Northern Ave., Boston, MA 42.34806 -71.03673

Esistono diversi modi per ottenere shapefile di codice postale. La più semplice è probabilmente le aree di tabulazione del codice postale dello US Census Bureau, che sono simili se non esattamente uguali ai confini del servizio postale degli Stati Uniti.

Puoi scaricare un file ZCTA direttamente dallo US Census Bureau, ma è un file per l'intero paese. Fallo solo se non ti dispiace un file di dati di grandi dimensioni. 

Un posto per scaricare un file ZCTA per un singolo stato è Census Reporter. Cerca i dati per stato, come la popolazione, quindi aggiungi il codice postale all'area geografica e scegli i dati di download come shapefile.

Potrei decomprimere manualmente il mio file scaricato, ma è più facile in R. Qui utilizzo la unzip()funzione di base R su un file scaricato e decomprimilo in una sottodirectory del progetto chiamata ma_zip_shapefile. Questo junkpaths = TRUEargomento dice che non voglio decomprimere l'aggiunta di un'altra sottodirectory basata sul nome del file zip.

unzip ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

sovrascrivi = TRUE)

Importazione e analisi geospaziale con sf

Ora finalmente un po 'di lavoro geospaziale. Importerò lo shapefile in R usando la st_read()funzione sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Lettura del layer `acs2017_5yr_B01003_86000US02648 'dalla sorgente dati` /Users/smachlis/Documents/MoreWipyr86/shape7_semplicemente con il driver di raccolta' ESBRI_Sshape048 ' caratteristiche e 4 campi # tipo di geometria: MULTIPOLYGON # dimensione: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs

Ho incluso la risposta della console durante l'esecuzione st_read()perché sono visualizzate alcune informazioni: epsg. Questo dice quale sistema di riferimento di coordinate è stato utilizzato per creare il file . Qui era 4326. Senza entrare troppo in profondità nelle erbacce, un epsg indica sostanzialmente  quale sistema è stato utilizzato per tradurre le aree di un globo tridimensionale - la Terra - in coordinate bidimensionali (latitudine e longitudine) . Questo è importante perché ci sono molti diversi sistemi di riferimento di coordinate. Voglio che i miei poligoni di codice postale e punti di indirizzo utilizzino lo stesso, in modo che si allineino correttamente.

Nota: questo file include un poligono per l'intero stato del Massachusetts, di cui non ho bisogno. Quindi filtrerò quella riga del Massachusetts con

zipcode_geo <- dplyr :: filter (zipcode_geo,

nome! = "Massachusetts")

Mappatura dello shapefile con tmap

La mappatura dei dati del poligono non è necessaria, ma è un bel controllo del mio shapefile per vedere se la geometria è quella che mi aspetto. Puoi fare un grafico veloce di un oggetto sf con la funzione di tmap qtm()(abbreviazione di quick theme map).

qtm (zipcode_geo) +

tm_legend (show = FALSE)

Schermate scattate da Sharon Machlis,

E sembra che io abbia effettivamente la geometria del Massachusetts con poligoni che potrebbero essere codici postali.

Successivamente voglio utilizzare i dati dell'indirizzo geocodificato. Questo è attualmente un semplice data frame, ma deve essere convertito in un oggetto geospaziale sf con il giusto sistema di coordinate.

Possiamo farlo con la st_as_sf()funzione sf . (Nota: le funzioni del pacchetto sf che operano su dati spaziali iniziano con st_, che sta per "spaziale" e "temporale".)

st_as_sf()richiede diversi argomenti. Nel codice seguente, il primo argomento è l'oggetto da trasformare: i miei indirizzi geocodificati. Il secondo vettore di argomento indica alla funzione quali colonne hanno i valori x (longitudine) e y (latitudine). Il terzo imposta il sistema di riferimento delle coordinate su 4326, quindi è lo stesso dei poligoni del codice postale.

point_geo <- st_as_sf (geocoded_addresses,

coords = c (x = "lon", y = "lat"),

crs = 4326)

Geospatial si unisce con sf

Ora che ho impostato i miei due set di dati, calcolare il codice postale per ogni indirizzo è facile con la st_join()funzione sf . La sintassi:

st_join (point_sf_object, polygon_sf_object, join = join_type)

In this example, I want to run st_join() on the geocoded points first and the ZIP code polygons second. It’s a so-called left join format: All points in the first data (geocoded addresses) are included, but only points in the second (ZIP code) data that match. Finally, my join type is st_within, since I want the match to be points within. 

my_results <- st_join(point_geo, zipcode_geo,

join = st_within)

That’s it! Now if I look at my results by printing out several of the most important columns, you”ll see each address has a ZIP code (in the “name” column). 

print(my_results[,c("query", "name", "geometry")])

# Semplice raccolta di caratteristiche con 2 caratteristiche e 2 campi # tipo di geometria: POINT # dimensione: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # nome query geometria # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Mappatura di punti e poligoni con tmap

Se desideri mappare i punti e i poligoni, ecco un modo per farlo con tmap:

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (my_results) +

tm_bubbles (col = "red", size = 0.25)

Schermata di Sharon Machlis,

Vuoi altri suggerimenti R? Vai alla pagina "Fai di più con R"!