Hur man gör rumsanalys i R med sf

Var röstar du? Vem är ni lagstiftare? Vad är ditt postnummer? Dessa frågor har något geografiskt gemensamt: Svaret innebär att bestämma vilken polygon en punkt faller inom.

Sådana beräkningar görs ofta med specialiserad GIS-programvara. Men de är också lätta att göra i R. Du behöver tre saker:

  1. Ett sätt att geokoda adresser för att hitta latitud och longitud; 
  2. Formfiler som beskriver polygongränser för postnummer; och 
  3. SF-paketet.

För geokodning använder jag vanligtvis geocod.io API. Det är gratis för 2500 uppslag om dagen och har ett trevligt R-paket, men du behöver en (gratis) API-nyckel för att använda den. För att komma runt den lite komplexiteten för den här artikeln använder jag det kostnadsfria Open Source Map Nominatim API. Det kräver ingen nyckel. Paketet tmaptools har en funktion,, geocode_OSM()att använda det API: et.

Importera och förbereda geospatial data

Jag använder paketen sf, tmaptools, tmap och dplyr. Om du vill följa med, ladda var och en med pacman::p_load()eller installera någon som ännu inte finns i ditt system med install.packages(), ladda sedan var och en med library().

För det här exemplet skapar jag en vektor med två adresser, vårt kontor i Framingham, Massachusetts, och RStudio-kontoret i Boston.

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

"250 Northern Ave., Boston, MA")

Geokodning är enkelt med geocode_OSM. Du kan se resultaten genom att skriva ut de tre första kolumnerna inklusive latitud och longitud:

geokodade_adresser <- geokod_OSM (adresser)

skriva ut (geokodade_adresser [, 1: 3])

fråga lat lon

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

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

Det finns flera sätt att hämta formfiler för postnummer. Det enklaste är förmodligen US Census Bureau: s postnummer Tabellområden, som liknar om inte exakt samma som US Postal Service-gränser.

Du kan ladda ner en ZCTA-fil direkt från US Census Bureau, men det är en fil för hela landet. Gör det bara om du inte har något emot en stor datafil. 

En plats att ladda ner en ZCTA-fil för en enda stat är Census Reporter. Sök efter data efter tillstånd, t.ex. befolkning, och lägg sedan till postnummer i geografin och välj nedladdningsdata som formfil.

Jag kunde packa upp min nedladdade fil manuellt, men det är lättare i R. Här använder jag bas R- unzip()funktionen på en nedladdad fil och packar upp den till en underkatalog projekt med namnet ma_zip_shapefile. Det junkpaths = TRUEargument säger att jag inte vill packa lägga till en annan underkatalog baserad på namnet på zip-filen.

packa upp ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", skräpvägar = SANT,

skriv över = SANT)

Geospatial import och analys med sf

Nu äntligen lite geospatialt arbete. Jag importerar shapefilen till R med sf: s st_read()funktion.

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Läslager `acs2017_5yr_B01003_86000US02648 'från datakällan` /Users/smachlis/Documents/MoreWithR/ma_zip_shop_sh funktioner och 4 fält # geometrityp: 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

Jag har inkluderat konsolsvaret när det körs st_read()eftersom det visas viss information där: epsg. Det säger vilket koordinatsreferenssystem som användes för att skapa filen . Här var det 4326. Utan att komma för djupt in i ogräset indikerar en epsg i princip  vilket system som användes för att översätta områden på en tredimensionell jordklot - jorden - till tvådimensionella koordinater (latitud och longitud) . Detta är viktigt eftersom det finns många olika referenssystem för koordinater. Jag vill att mina polygoner och adresspunkter ska använda samma, så att de stämmer ordentligt.

Obs! Den här filen innehåller en polygon för hela Massachusetts, vilket jag inte behöver. Så jag ska filtrera bort den Massachusetts-raden med

zipcode_geo <- dplyr :: filter (zipcode_geo,

namn! = "Massachusetts")

Kartlägga formfilen med tmap

Det är inte nödvändigt att kartlägga polygondata, men det är en trevlig kontroll av min shapefile för att se om geometrin är vad jag förväntar mig. Du kan göra en snabb plottning av ett sf-objekt med tmaps qtm()(kort för snabb temakarta) funktion.

qtm (postnummer_geo) +

tm_legend (visa = FALSE)

Skärmar tagna av Sharon Machlis,

Och det ser ut som om jag verkligen har Massachusetts-geometri med polygoner som kan vara postnummer.

Därefter vill jag använda geokodad adressinformation. Detta är för närvarande en vanlig dataram, men den måste konverteras till ett sf geospatialt objekt med rätt koordinatsystem.

Vi kan göra det med sf: s st_as_sf()funktion. (Obs: sf-paketfunktioner som fungerar på rumsliga data börjar med st_, vilket står för "rumslig" och "tidsmässig.")

st_as_sf()tar flera argument. I koden nedan är det första argumentet objektet att transformera - mina geokodade adresser. Den andra argumentvektorn berättar funktionen vilka kolumner som har värdena x (longitud) och y (latitud). Den tredje ställer in koordinatsreferenssystemet till 4326, så det är detsamma som mina postnummer polygoner.

point_geo <- st_as_sf (geokodade adresser,

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

crs = 4326)

Geospatial ansluter sig till sf

Nu när jag har ställt in mina två datamängder är det enkelt att beräkna postnummer för varje adress med sf: s st_join()funktion. Syntaxen:

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

# Enkel funktionssamling med 2 funktioner och 2 fält # geometrityp: PUNKT # 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 # frågens namngeometri # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Kartläggningspunkter och polygoner med tmap

Om du vill mappa punkter och polygoner, här är ett sätt att göra det med tmap:

tm_shape (postnummer_geo) +

tm_fill () +

tm_shape (min_resultat) +

tm_bubbles (kol = "röd", storlek = 0,25)

Skärmdump av Sharon Machlis,

Vill du ha fler R-tips? Gå till sidan "Gör mer med R"!