Wie man eine räumliche Analyse in R mit sf durchführt

Wo stimmen Sie ab? Wer sind Sie Gesetzgeber? Wie ist Ihre Postleitzahl? Diese Fragen haben räumlich etwas gemeinsam: Die Antwort besteht darin, zu bestimmen, in welches Polygon ein Punkt fällt.

Solche Berechnungen werden häufig mit einer speziellen GIS-Software durchgeführt. Aber sie sind auch in R einfach zu machen. Sie brauchen drei Dinge:

  1. Eine Möglichkeit, Adressen zu geocodieren, um Längen- und Breitengrade zu ermitteln. 
  2. Shapefiles, die Postleitzahl-Polygongrenzen umreißen; und 
  3. Das sf-Paket.

Für die Geokodierung verwende ich normalerweise die geocod.io-API. Es ist kostenlos für 2.500 Suchvorgänge pro Tag und hat ein schönes R-Paket, aber Sie benötigen einen (kostenlosen) API-Schlüssel, um es zu verwenden. Um diese Komplexität für diesen Artikel zu umgehen, verwende ich die kostenlose Open-Source-Open Street Map Nominatim-API. Es ist kein Schlüssel erforderlich. Das tmaptools-Paket hat die Funktion, geocode_OSM()diese API zu verwenden.

Geodaten importieren und vorbereiten

Ich werde die Pakete sf, tmaptools, tmap und dplyr verwenden. Wenn Sie mitmachen möchten, laden Sie jede mit pacman::p_load()oder installieren Sie noch nicht auf Ihrem System mit install.packages(), und laden Sie dann jede mit library().

In diesem Beispiel erstelle ich einen Vektor mit zwei Adressen, unserem Büro in Framingham, Massachusetts, und dem RStudio-Büro in Boston.

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

"250 Northern Ave., Boston, MA")

Die Geokodierung ist mit geocode_OSM unkompliziert. Sie können die Ergebnisse anzeigen, indem Sie die ersten drei Spalten einschließlich Breiten- und Längengrad ausdrucken:

geocodierte_Adressen <- geocode_OSM (Adressen)

print (geocoded_addresses [, 1: 3])

Abfrage lat lon

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

# 2 250 Northern Ave., Boston, MA 42.34806 -71.03673

Es gibt verschiedene Möglichkeiten, Postleitzahl-Shapefiles abzurufen. Am einfachsten sind wahrscheinlich die Tabellierungsbereiche der Postleitzahl des US Census Bureau, die den Grenzen des US-Postdienstes ähneln, wenn nicht sogar genau diesen entsprechen.

Sie können eine ZCTA-Datei direkt vom US Census Bureau herunterladen, sie ist jedoch eine Datei für das gesamte Land. Tun Sie das nur, wenn Ihnen eine große Datendatei nichts ausmacht. 

Ein Ort zum Herunterladen einer ZCTA-Datei für einen einzelnen Status ist Census Reporter. Suchen Sie nach Daten nach Bundesland, z. B. Bevölkerung, und fügen Sie der Geografie eine Postleitzahl hinzu. Wählen Sie dann Daten als Shapefile herunterladen.

Ich könnte meine heruntergeladene Datei manuell entpacken, aber in R ist es einfacher. Hier verwende ich die unzip()Funktion von base R für eine heruntergeladene Datei und entpacke sie in ein Projekt-Unterverzeichnis namens ma_zip_shapefile. Dieses junkpaths = TRUEArgument besagt, dass ich nicht entpacken möchte, indem ich ein weiteres Unterverzeichnis basierend auf dem Namen der Zip-Datei hinzufüge.

entpacken ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

überschreiben = WAHR)

Geodatenimport und Analyse mit sf

Nun endlich etwas Geodatenarbeit. Ich werde das Shapefile mit der st_read()Funktion von sf in R importieren .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Leseebene "acs2017_5yr_B01003_86000US02648" aus der Datenquelle "/Users/smachlis/Documents/MoreWithR/ma_p" Features und 4 Felder # Geometrietyp: MULTIPOLYGON # Dimension: XY # bbox: xmin: -73,50821 ymin: 41,18705 xmax: -69,85886 ymax: 42,95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs

Ich habe die Konsolenantwort beim Ausführen st_read()eingefügt , da dort einige Informationen angezeigt werden: das epsg. Das heißt, welches Koordinatenreferenzsystem zum Erstellen der Datei verwendet wurde . Hier war es 4326. Ohne zu tief in das Unkraut einzudringen, gibt ein epsg im Grunde an, mit  welchem ​​System Gebiete auf einem dreidimensionalen Globus - der Erde - in zweidimensionale Koordinaten (Breiten- und Längengrade) übersetzt wurden . Dies ist wichtig, da es viele verschiedene Koordinatenreferenzsysteme gibt. Ich möchte, dass meine Postleitzahl-Polygone und Adresspunkte dasselbe verwenden, damit sie richtig ausgerichtet sind.

Hinweis: Diese Datei enthält zufällig ein Polygon für den gesamten Bundesstaat Massachusetts, das ich nicht benötige. Also werde ich diese Massachusetts-Reihe mit herausfiltern

Postleitzahl_geo <- dplyr :: filter (Postleitzahl_geo,

name! = "Massachusetts")

Mapping des Shapefiles mit tmap

Das Zuordnen der Polygondaten ist nicht erforderlich, aber es ist eine schöne Überprüfung meines Shapefiles, um festzustellen, ob die Geometrie meinen Erwartungen entspricht. Mit der Funktion tmap qtm()(kurz für Quick Theme Map) können Sie ein SF-Objekt schnell zeichnen .

qtm (Postleitzahl_geo) +

tm_legend (show = FALSE)

Bildschirme von Sharon Machlis,

Und es sieht so aus, als hätte ich tatsächlich eine Massachusetts-Geometrie mit Polygonen, die Postleitzahlen sein könnten.

Als nächstes möchte ich die geokodierten Adressdaten verwenden. Dies ist derzeit ein einfacher Datenrahmen, der jedoch in ein sf-Geodatenobjekt mit dem richtigen Koordinatensystem konvertiert werden muss.

Wir können das mit der st_as_sf()Funktion von sf machen . (Hinweis: sf-Paketfunktionen, die mit Geodaten arbeiten, beginnen mit st_, was für "räumlich" und "zeitlich" steht.)

st_as_sf()nimmt mehrere Argumente. Im folgenden Code ist das erste Argument das zu transformierende Objekt - meine geokodierten Adressen. Der zweite Argumentvektor teilt der Funktion mit, welche Spalten die Werte x (Länge) und y (Breite) haben. Das dritte setzt das Koordinatenreferenzsystem auf 4326, es ist also dasselbe wie meine Postleitzahl-Polygone.

point_geo <- st_as_sf (geocodierte_Adressen,

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

crs = 4326)

Geospatial schließt sich sf an

Nachdem ich meine beiden Datensätze eingerichtet habe, ist die Berechnung der Postleitzahl für jede Adresse mit der st_join()Funktion von sf ganz einfach . Die Syntax:

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

# Einfache Feature-Sammlung mit 2 Features und 2 Feldern # Geometrietyp: POINT # Dimension: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # Abfrage Name Geometrie # 1 492 Alter Connecticut-Pfad, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2 250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Zuordnungspunkte und Polygone mit tmap

Wenn Sie die Punkte und Polygone zuordnen möchten, haben Sie folgende Möglichkeiten, dies mit tmap zu tun:

tm_shape (Postleitzahl_geo) +

tm_fill () +

tm_shape (my_results) +

tm_bubbles (col = "rot", Größe = 0,25)

Screenshot von Sharon Machlis,

Willst du mehr R-Tipps? Gehen Sie zur Seite "Mehr mit R machen"!