Hvor stemmer du? Hvem er dere lovgivere? Hva er postnummeret ditt? Disse spørsmålene har noe geografisk til felles: Svaret innebærer å bestemme hvilken polygon et punkt faller innenfor.
Slike beregninger gjøres ofte med spesialisert GIS-programvare. Men de er også enkle å gjøre i R. Du trenger tre ting:
- En måte å geokode adresser for å finne breddegrad og lengdegrad;
- Formfil som skisserer polygongrenser for postnummer; og
- SF-pakken.
For geokoding bruker jeg vanligvis geocod.io API. Det er gratis for 2500 oppslag om dagen og har en fin R-pakke, men du trenger en (gratis) API-nøkkel for å bruke den. For å omgå den litt kompleksiteten for denne artikkelen, bruker jeg gratis, åpen kildekode Open Street Map Nominatim API. Det krever ikke en nøkkel. Tmaptools-pakken har en funksjon, geocode_OSM (), for å bruke API-en.
Import og prepping geospatial data
Jeg bruker pakkene sf, tmaptools, tmap og dplyr. Hvis du vil følge med, laster du hver med pacman :: p_load () eller installer noen som ennå ikke er på systemet ditt med install.packages (), last deretter hver med bibliotek().
For dette eksemplet lager jeg en vektor med to adresser, vårt kontor i Framingham, Massachusetts, og RStudio-kontoret i Boston.
adresser <- c ("492 Old Connecticut Path, Framingham, MA","250 Northern Ave., Boston, MA")
Geokoding er grei med geocode_OSM. Du kan se resultatene ved å skrive ut de tre første kolonnene, inkludert breddegrad og lengdegrad:
geokodede adresser <- geokode_OSM (adresser)skriv ut (geokodede_adresser [, 1: 3])
spørring lat lon
# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105
# 2 250 Northern Ave., Boston, MA 42.34806 -71.03673
Det er flere måter å få postnummerformer på. Det enkleste er sannsynligvis US Census Bureau's ZIP Code Tabulation Areas, som ligner på om ikke akkurat det samme som U.S. Postal Service-grenser.
Du kan laste ned en ZCTA-fil direkte fra US Census Bureau, men det er en fil for hele landet. Gjør det bare hvis du ikke har noe imot en stor datafil.
Et sted å laste ned en ZCTA-fil for en enkelt stat er Census Reporter. Søk etter data etter tilstand, for eksempel populasjon, og legg deretter til postnummer i geografien, og velg nedlastingsdata som formfil.
Jeg kunne pakke ut den nedlastede filen min manuelt, men det er lettere i R. Her bruker jeg base R-er pakke ut () funksjon på en nedlastet fil, og pakke den ut til en prosjektunderkatalog med navnet ma_zip_shapefile. At søppelstier = SANT argumentet sier at jeg ikke vil pakke ut å legge til en annen underkatalog basert på navnet på zip-filen.
pakke ut ("data / acs2017_5yr_B01003_86000US02648.zip",exdir = "ma_zip_shapefile", søppelstier = SANT,
overskriv = SANT)
Geospatial import og analyse med sf
Nå endelig noe geospatialt arbeid. Jeg importerer shapefilen til R ved hjelp av sf st_read () funksjon.
zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Leselag `acs2017_5yr_B01003_86000US02648 'fra datakilde' /Users/smachlis/Documents/MoreWithR/ma_zip_shop_shop funksjoner og 4 felt # geometritype: MULTIPOLYGON # dimensjon: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defsJeg har tatt med konsollresponsen når du kjører st_read () fordi det vises litt informasjon der: epsg. Det sier hvilket koordinatsystem som ble brukt til å lage filen. Her var det 4326. Uten å komme for dypt inn i ugresset, indikerer en epsg i utgangspunktethvilket system ble brukt til å oversette områder på en tredimensjonal klode - jorden - til todimensjonale koordinater (breddegrad og lengdegrad). Dette er viktig fordi det er en mye av forskjellige koordinatsystemer. Jeg vil at postnummer polygoner og adressepunkter skal bruke den samme, slik at de stiller opp riktig.
Merk: Denne filen inneholder tilfeldigvis en polygon for hele staten Massachusetts, som jeg ikke trenger. Så jeg filtrerer ut den Massachusetts-raden med
zipcode_geo <- dplyr :: filter (zipcode_geo,navn! = "Massachusetts")
Kartlegge shapefilen med tmap
Kartlegging av polygondataene er ikke nødvendig, men det er en fin sjekk av shapefilen min for å se om geometrien er det jeg forventer. Du kan gjøre et raskt plott av et sf-objekt med tmap’s qtm () (kort for rask temakart) -funksjon.
qtm (postnummer_geo) +tm_legend (show = FALSE)
Skjermbilder tatt av Sharon Machlis, Og det ser ut til at jeg faktisk har Massachusetts-geometri med polygoner som kan være postnummer.
Deretter vil jeg bruke geokodede adressedata. Dette er for tiden en vanlig dataramme, men den må konverteres til et geospatialt objekt med riktig koordinatsystem.
Vi kan gjøre det med sf’s st_as_sf () funksjon. (Merk: sf-pakkefunksjoner som fungerer på geodata starter med st_, som står for "romlig" og "tidsmessig.")
st_as_sf () tar flere argumenter. I koden nedenfor er det første argumentet objektet som skal transformeres - de geokodede adressene mine. Den andre argumentvektoren forteller funksjonen hvilke kolonner som har x (lengdegrad) og y (breddegrad). Den tredje setter koordinatsystemet til 4326, så det er det samme som postnummerpolygonene mine.
point_geo <- st_as_sf (geokodede adresser,koordiner = c (x = "lon", y = "lat"),
crs = 4326)
Geospatial slutter seg til sf
Nå som jeg har satt opp mine to datasett, er det enkelt å beregne postnummeret for hver adresse med sf-er st_join () funksjon. Syntaksen:
st_join (point_sf_object, polygon_sf_object, join = join_type)
I dette eksemplet vil jeg løpe st_join () på geokodede punkter først og postnummer polygoner andre. Det er et såkalt left join format: Alle poeng i de første dataene (geokodede adresser) er inkludert, men bare poeng i de andre (postnummer) dataene som samsvarer. Endelig er min type st_intern, siden jeg vil at kampen skal være poeng innenfor.
mine_resultater <- st_join (point_geo, zipcode_geo,bli med = st_within)
Det er det! Nå hvis jeg ser på resultatene mine ved å skrive ut flere av de viktigste kolonnene, vil du se at hver adresse har et postnummer (i "navn" -kolonnen).
skriv ut (mine_resultater [, c ("spørring", "navn", "geometri")])# Enkel funksjonssamling med 2 funksjoner og 2 felt # geometritype: PUNKT # dimensjon: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # query name geometry # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)
Kartleggingspunkter og polygoner med tmap
Hvis du vil kartlegge punktene og polygonene, kan du gjøre det med tmap:
tm_shape (postnummer_geo) +tm_fill () +
tm_shape (mine_resultater) +
tm_bubbles (col = "rød", størrelse = 0,25)
Skjermbilde av Sharon Machlis, Vil du ha flere R-tips? Gå til siden "Gjør mer med R"!
