Hvor stemmer du? Hvem er I lovgivere? Hvad er dit postnummer? Disse spørgsmål har noget geospatielt til fælles: Svaret indebærer at bestemme, hvilken polygon et punkt falder inden for.
Sådanne beregninger udføres ofte med specialiseret GIS-software. Men de er også nemme at gøre i R. Du har brug for tre ting:
- En måde at geokode adresser for at finde bredde og længdegrad;
- Shapefiles, der skitserer postnummer polygon grænser; og
- SF-pakken.
Til geokodning bruger jeg normalt geocod.io API. Det er gratis til 2.500 opslag om dagen og har en dejlig R-pakke, men du har brug for en (gratis) API-nøgle for at bruge den. For at omgå den smule af kompleksitet i denne artikel bruger jeg den gratis open source Open Street Map Nominatim API. Det kræver ikke en nøgle. Tmaptools-pakken har en funktion, geocode_OSM (), for at bruge denne API.
Import og prepping geospatiale data
Jeg bruger pakkerne sf, tmaptools, tmap og dplyr. Hvis du vil følge med, skal du indlæse hver med pacman :: p_load () eller installer nogen, der endnu ikke er på dit system med install.packages (), ilæg derefter hver med bibliotek().
Til dette eksempel opretter jeg en vektor med to adresser, vores kontor i Framingham, Massachusetts, og RStudio-kontoret i Boston.
adresser <- c ("492 Old Connecticut Path, Framingham, MA","250 Northern Ave., Boston, MA")
Geokodning er ligetil med geocode_OSM. Du kan se resultaterne ved at udskrive de første tre kolonner inklusive breddegrad og længdegrad:
geokodede adresser <- geokode_OSM (adresser)udskriv (geokodede_adresser [, 1: 3])
forespørgsel lat lon
# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105
# 2250 Northern Ave., Boston, MA 42.34806 -71.03673
Der er flere måder at få postnummerformularer på. Det nemmeste er sandsynligvis det amerikanske folketællingsbureaus postnummeroversigtsområder, der ligner, hvis ikke nøjagtigt det samme som U.S. Postal Service-grænser.
Du kan downloade en ZCTA-fil direkte fra U.S.Census Bureau, men det er en fil for hele landet. Gør det kun, hvis du ikke har noget imod en stor datafil.
Et sted at downloade en ZCTA-fil til en enkelt stat er Census Reporter. Søg efter data efter tilstand, f.eks. Population, og tilføj derefter postnummer til geografien, og vælg download-data som shapefile.
Jeg kunne udpakke min downloadede fil manuelt, men det er lettere i R. Her bruger jeg base R'er pakke ud () funktion på en downloadet fil og pakke den ud til en projektunderkatalog med navnet ma_zip_shapefile. At junkpaths = SAND argument siger, at jeg ikke vil pakke ud tilføjelse af en anden underkatalog baseret på navnet på zip-filen.
pakke ud ("data / acs2017_5yr_B01003_86000US02648.zip",exdir = "ma_zip_shapefile", junkpaths = SAND,
overskriv = SAND)
Geospatial import og analyse med sf
Nu endelig noget geospatialt arbejde. Jeg importerer shapefilen til R ved hjælp af sf'er st_read () fungere.
zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Læselag `acs2017_5yr_B01003_86000US02648 'fra datakilde' /Users/smachlis/Documents/MoreWithR/ma_zip_shop_shop_shop funktioner og 4 felter # geometri type: 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_defsJeg har medtaget konsolresponset, når det kører st_read () fordi der vises nogle oplysninger der: epsg. Det siger hvilket koordinatsystem der blev brugt til at oprette filen. Her var det 4326. Uden at komme for dybt ind i ukrudtet indikerer en epsg grundlæggendehvilket system blev brugt til at oversætte områder på en tredimensionel klode - Jorden - til todimensionale koordinater (breddegrad og længdegrad). Dette er vigtigt, fordi der er en masse af forskellige koordinatsystemer. Jeg vil have, at mine postnummer polygoner og adressepunkter skal bruge den samme, så de stemmer overens.
Bemærk: Denne fil indeholder tilfældigvis en polygon for hele staten Massachusetts, som jeg ikke har brug for. Så jeg filtrerer den Massachusetts-række ud med
postnummer_geo <- dplyr :: filter (postnummer_geo,navn! = "Massachusetts")
Kortlægning af shapefilen med tmap
Kortlægning af polygondata er ikke nødvendigt, men det er en god kontrol af min shapefile for at se, om geometrien er, hvad jeg forventer. Du kan lave et hurtigt plot af et sf-objekt med tmap's qtm () (forkortelse for hurtig temakort) -funktion.
qtm (postnummer_geo) +tm_legend (show = FALSE)
Skærme skudt af Sharon Machlis, Og det ser ud til, at jeg faktisk har Massachusetts-geometri med polygoner, der kunne være postnumre.
Dernæst vil jeg bruge de geokodede adressedata. Dette er i øjeblikket en almindelig dataramme, men den skal konverteres til et sf geospatialt objekt med det rigtige koordinatsystem.
Vi kan gøre det med sf'er st_as_sf () fungere. (Bemærk: sf-pakkefunktioner, der fungerer på geodata, starter med st_, som står for "rumlig" og "tidsmæssig.")
st_as_sf () tager flere argumenter. I koden nedenfor er det første argument det objekt, der skal transformeres - mine geokodede adresser. Den anden argumentvektor fortæller funktionen, hvilke kolonner der har værdierne x (længdegrad) og y (breddegrad). Den tredje sætter koordinatsystemet til 4326, så det er det samme som mine postnummer polygoner.
point_geo <- st_as_sf (geokodede_adresser,koordiner = c (x = "lon", y = "lat"),
crs = 4326)
Geospatial slutter sig til sf
Nu hvor jeg har oprettet mine to datasæt, er det nemt at beregne postnummeret for hver adresse med sf'er st_join () fungere. Syntaksen:
st_join (point_sf_object, polygon_sf_object, join = join_type)
I dette eksempel vil jeg køre st_join () på de geokodede punkter først og postnummer polygoner andet. Det er et såkaldt left join format: Alle punkter i de første data (geokodede adresser) er inkluderet, men kun punkter i de andet (postnummer) data, der matcher. Endelig er min tilslutningstype st_intern, da jeg ønsker, at kampen skal være point indeni.
mine_resultater <- st_join (point_geo, postnummer_geo,slutte sig til = st_within)
Det er det! Hvis jeg nu ser på mine resultater ved at udskrive flere af de vigtigste kolonner, vil du se, at hver adresse har et postnummer (i kolonnen "navn").
print (mine_resultater [, c ("forespørgsel", "navn", "geometri")])# Enkel funktionssamling med 2 funktioner og 2 felter # geometri type: 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 # forespørgselsnavn geometri # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)
Kortlægningspunkter og polygoner med tmap
Hvis du vil kortlægge punkterne og polygonerne, er der en måde at gøre det på med tmap:
tm_shape (postnummer_geo) +tm_fyld () +
tm_shape (mine_resultater) +
tm_bubbles (col = "rød", størrelse = 0,25)
Screen shot af Sharon Machlis, Vil du have flere R-tip? Gå til siden "Gør mere med R"!
