Kde hlasujete? Kdo jste zákonodárci? Jaké je vaše PSČ? Tyto otázky mají něco geoprostorového společného: Odpověď zahrnuje určení toho, do kterého polygonu bod spadá.
Takové výpočty se často provádějí pomocí specializovaného softwaru GIS. Ale také se snadno dělají v R. Potřebujete tři věci:
- Způsob geokódování adres pro zjištění zeměpisné šířky a délky;
- Shapefiles, které vymezují hranice polygonů PSČ; a
- Balíček sf.
Pro geokódování obvykle používám API geocod.io. Je zdarma pro 2 500 vyhledávání denně a má pěkný balíček R, ale k jeho použití potřebujete (bezplatný) klíč API. Abych pro tento článek obešel tuto složitost, použiji bezplatné open-source rozhraní Open Street Map Nominatim API. To nevyžaduje klíč. Balíček tmaptools má funkci, geocode_OSM (), používat toto API.
Import a příprava geoprostorových dat
Budu používat balíčky sf, tmaptools, tmap a dplyr. Pokud chcete sledovat, načtěte každý s pacman :: p_load () nebo nainstalujte všechny, které dosud ve vašem systému nejsou install.packages (), poté načtěte každý knihovna().
V tomto příkladu vytvořím vektor se dvěma adresami, naší kancelář ve Framinghamu v Massachusetts a kancelář RStudio v Bostonu.
adresy <- c ("492 Old Connecticut Path, Framingham, MA",„250 Northern Ave., Boston, MA“)
Geokódování je s geocode_OSM jednoduché. Výsledky můžete zobrazit vytištěním prvních tří sloupců včetně zeměpisné šířky a délky:
geocoded_addresses <- geocode_OSM (adresy)tisk (geokódované_adresy [, 1: 3])
dotaz lat lon
1 492 Old Connecticut Path, Framingham, MA 42.31348 - 71.39105
# 2 250 Northern Ave., Boston, MA 42,34806 - 71,03673
Existuje několik způsobů, jak získat tvarové soubory PSČ. Nejjednodušší je pravděpodobně tabelační oblast poštovního směrovacího čísla amerického sčítání lidu, která je podobná, ne-li úplně stejná, jako hranice americké poštovní služby.
Soubor ZCTA si můžete stáhnout přímo z amerického sčítání lidu, ale je to soubor pro celou zemi. Udělejte to, pouze pokud vám nevadí velký datový soubor.
Jedno místo pro stažení souboru ZCTA pro jeden stát je Census Reporter. Vyhledejte libovolná data podle stavu, například počet obyvatel, a poté přidejte PSČ do geografie a zvolte stahovaná data jako obrazec.
Stažený soubor jsem mohl rozbalit ručně, ale v R. je to jednodušší. Tady používám základní R. rozbalit () funkce na staženém souboru a rozbalte jej do podadresáře projektu s názvem ma_zip_shapefile. Že junkpaths = PRAVDA Argument říká, že nechci rozbalit přidání dalšího podadresáře na základě názvu souboru zip.
unzip ("data / acs2017_5yr_B01003_86000US02648.zip",exdir = "ma_zip_shapefile", junkpaths = PRAVDA,
přepsat = PRAVDA)
Geoprostorový import a analýza se sf
Nyní konečně geoprostorové práce. Importuji tvarový soubor do R pomocí sf st_read () funkce.
zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # čtecí vrstva `acs2017_5yr_B01003_86000US02648 'ze zdroje dat' /Users/smachlis/Documents/M_W_Ro__8_r_ ' vlastnosti a 4 pole # typ geometrie: MULTIPOLYGON # rozměr: XY # bbox: xmin: -73,50821 ymin: 41,18705 xmax: -69,85886 ymax: 42,95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defsPři spuštění jsem zahrnul odpověď konzoly st_read () protože se tam zobrazují některé informace: epsg. To říká jaký souřadnicový referenční systém byl použit k vytvoření souboru. Tady to bylo 4326. Aniž bychom se dostali příliš hluboko do plevelů, epsg v podstatě naznačujejaký systém byl použit k překladu oblastí na trojrozměrném světě - Zemi - na dvourozměrné souřadnice (zeměpisná šířka a délka). To je důležité, protože existují hodně různých souřadnicových referenčních systémů. Chci, aby mé polygony a adresy v PSČ používaly stejný, aby se správně seřadily.
Poznámka: Tento soubor obsahuje polygon pro celý stát Massachusetts, což nepotřebuji. Odfiltruji tedy ten řádek v Massachusetts
zipcode_geo <- dplyr :: filter (zipcode_geo,name! = "Massachusetts")
Mapování obrazce pomocí tmap
Mapování polygonových dat není nutné, ale je to pěkná kontrola mého obrazce, zda geometrie odpovídá očekávání. Pomocí tmap můžete rychle vykreslit objekt sf qtm () (zkratka pro rychlou mapu témat).
qtm (zipcode_geo) +tm_legend (show = FALSE)
Obrazovky natočené Sharon Machlis, A vypadá to, že opravdu mám Massachusettskou geometrii s polygony, které by mohly být PSČ.
Dále chci použít data geokódované adresy. Toto je v současné době prostý datový rámec, ale je třeba jej převést na geoprostorový objekt s vhodným souřadným systémem.
Můžeme to udělat se sf st_as_sf () funkce. (Poznámka: Funkce balíčku sf, které fungují na prostorových datech, začínají Svatý_, což znamená „prostorové“ a „časové“.)
st_as_sf () trvá několik argumentů. V níže uvedeném kódu je prvním argumentem objekt, který se má transformovat - moje geokódované adresy. Druhý vektor argumentů říká funkci, které sloupce mají hodnoty x (zeměpisná délka) a y (zeměpisná šířka). Třetí nastavuje souřadnicový referenční systém na 4326, takže je stejný jako moje polygony s PSČ.
point_geo <- st_as_sf (geokódované_adresy,coords = c (x = "lon", y = "lat"),
crs = 4326)
Geoprostorové spojení se sf
Nyní, když jsem nastavil své dvě datové sady, je výpočet PSČ pro každou adresu pomocí sf snadný st_join () funkce. Syntaxe:
st_join (point_sf_object, polygon_sf_object, join = join_type)
V tomto příkladu chci spustit st_join () nejprve na geokódované body a na druhém polygony PSČ. Jedná se o takzvaný levý formát připojení: Všechno jsou zahrnuty body v prvních datech (geokódované adresy), ale pouze body v druhém (PSČ) údaje, které se shodují. Nakonec je můj typ spojení st_within, protože chci, aby zápas byl uvnitř.
my_results <- st_join (point_geo, zipcode_geo,join = st_within)
A je to! Nyní, když se podívám na své výsledky vytištěním několika nejdůležitějších sloupců, uvidíte, že každá adresa má PSČ (ve sloupci „název“).
print (my_results [, c ("dotaz", "jméno", "geometrie")])# Jednoduchá kolekce funkcí se 2 funkcemi a 2 poli # typ geometrie: POINT # rozměr: XY # bbox: xmin: -71,39105 ymin: 42,31348 xmax: -71,03673 ymax: 42,34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # geometrie názvu dotazu # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71,39105 42,31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71,03673 42,34806)
Mapování bodů a polygonů pomocí tmap
Pokud byste chtěli zmapovat body a polygony, je tu jeden způsob, jak to udělat pomocí tmap:
tm_shape (zipcode_geo) +tm_fill () +
tm_shape (my_results) +
tm_bubbles (col = "red", velikost = 0,25)
Snímek obrazovky Sharon Machlis, Chcete více R tipů? Přejděte na stránku „Udělejte více s R“!
