Kde volíte? Kto ste zákonodarcovia? Aké je vaše PSČ? Tieto otázky majú niečo geopriestorovo spoločné: Odpoveď spočíva v určení, do ktorého mnohouholníka bod spadá.
Takéto výpočty sa často vykonávajú so špecializovaným softvérom GIS. Ľahko sa však dajú robiť aj v R. Potrebujete tri veci:
- Spôsob geokódovania adries na zistenie zemepisnej šírky a dĺžky;
- Súbory tvaru, ktoré ohraničujú hranice polygónov PSČ; a
- Balíček sf.
Na geokódovanie zvyčajne používam rozhranie geocod.io API. Je zadarmo pre 2 500 vyhľadávaní denne a má pekný balík R. Na jeho použitie však potrebujete (bezplatný) kľúč API. Aby som tento článok obišiel tak trochu zložito, použijem bezplatné open source API Open Street Map Nominatim API. Nevyžaduje kľúč. Balík tmaptools má funkciu, geocode_OSM (), používať toto API.
Import a príprava geopriestorových údajov
Budem používať balíčky sf, tmaptools, tmap a dplyr. Ak chcete sledovať, každý načítajte pacman :: p_load () alebo nainštalujte všetky, ktoré ešte vo vašom systéme nie sú, install.packages (), potom načítajte každý s knižnica ().
Pre tento príklad vytvorím vektor s dvoma adresami, našu kanceláriu vo Framinghame v štáte Massachusetts a kanceláriu RStudio v Bostone.
adresy <- c ("492 Old Connecticut Path, Framingham, MA",„250 Northern Ave., Boston, MA“)
Geokódovanie je s geocode_OSM jednoduché. Výsledky si môžete pozrieť vytlačením prvých troch stĺpcov vrátane zemepisnej šírky a dĺžky:
geocoded_addresses <- geocode_OSM (adresy)tlač (geocoded_addresses [, 1: 3])
dopyt lat lon
1 492 Old Connecticut Path, Framingham, MA 42.31348 - 71.39105
# 2 250 Northern Ave., Boston, MA 42.34806 -71,03673
Existuje niekoľko spôsobov, ako získať súbory s tvarom PSČ. Najjednoduchšie sú pravdepodobne tabulačné oblasti poštového smerovacieho čísla amerického úradu pre sčítanie ľudu, ktoré sú podobné, ak nie úplne rovnaké, ako hranice americkej poštovej služby.
Môžete si stiahnuť súbor ZCTA priamo z Úradu pre sčítanie ľudu USA, ale je to súbor pre celú krajinu. Urobte to iba v prípade, že vám nevadí veľký dátový súbor.
Jedným z miest na stiahnutie súboru ZCTA pre jeden štát je Census Reporter. Vyhľadajte ľubovoľné údaje podľa stavu, napríklad počet obyvateľov, potom pridajte do geografického smerovacieho čísla PSČ a zvoľte sťahovacie údaje ako obrazcový súbor.
Stiahnutý súbor som mohol rozbaliť manuálne, ale v R. Je to jednoduchšie. Tu používam základné R. rozbaliť () funkciu na stiahnutom súbore a rozbaľte ho do podadresára projektu s názvom ma_zip_shapefile. To junkpaths = TRUE argument hovorí, že nechcem rozbaliť pridanie ďalšieho podadresára na základe názvu zip súboru.
unzip ("data / acs2017_5yr_B01003_86000US02648.zip",exdir = "ma_zip_shapefile", junkpaths = TRUE,
prepísať = TRUE)
Geopriestorový import a analýza so sf
Teraz konečne nejaké geopriestorové práce. Importujem tvarový súbor do R pomocou sf st_read () funkcia.
zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # vrstva na čítanie `acs2017_5yr_B01003_86000US02648 'zo zdroja údajov' /Users/smachlis/Documents/MoreWithR/ma_z_8_r_sh_r8 ' funkcie a 4 polia # typ geometrie: MULTIPOLYGON # rozmer: XY # bbox: xmin: -73,50821 ymin: 41,18705 xmax: -69,85886 ymax: 42,95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defsPri spustení som zahrnul odpoveď konzoly st_read () pretože sú tam zobrazené niektoré informácie: epsg. To hovorí aký súradnicový referenčný systém bol použitý na vytvorenie súboru. Tu to bolo 4326. Bez toho, aby sme sa dostali príliš hlboko do buriny, epsg v podstate naznačujeaký systém bol použitý na preklad oblastí na trojrozmernej planéte - Zemi - na dvojrozmerné súradnice (zemepisná šírka a dĺžka). To je dôležité, pretože existujú a veľa rôznych referenčných systémov súradníc. Chcem, aby moje polygóny s PSČ a adresné body používali ten istý, aby sa správne zoradili.
Poznámka: Tento súbor obsahuje polygón pre celý štát Massachusetts, čo nepotrebujem. Takže odfiltrujem ten rad z Massachusetts
zipcode_geo <- dplyr :: filter (zipcode_geo,name! = "Massachusetts")
Mapovanie tvarového súboru s tmap
Mapovanie údajov o mnohouholníku nie je potrebné, ale je to pekná kontrola môjho súboru tvaru, aby som zistil, či je geometria taká, akú očakávam. Smm môžete urobiť rýchle vykreslenie sf objektu qtm () (skratka pre rýchlu mapu tém).
qtm (zipcode_geo) +tm_legend (show = FALSE)
Obrazovky natočené Sharon Machlis, A vyzerá to, že skutočne mám Massachusettskú geometriu s mnohouholníkmi, ktoré by mohli byť PSČ.
Ďalej chcem použiť geokódované údaje o adrese. Toto je v súčasnosti obyčajný dátový rámec, ale je potrebné ho previesť na geopriestorový objekt s vhodným súradnicovým systémom.
Môžeme to urobiť pomocou sf’s st_as_sf () funkcia. (Poznámka: funkcie balíka sf, ktoré fungujú na priestorových dátach, začínajú na st_, čo znamená „priestorové“ a „časové“.)
st_as_sf () trvá niekoľko argumentov. V kóde nižšie je prvým argumentom objekt, ktorý sa má transformovať - moje geokódované adresy. Druhý vektor argumentov hovorí funkcii, ktoré stĺpce majú hodnoty x (zemepisná dĺžka) a y (zemepisná šírka). Tretia nastavuje referenčný súradnicový systém na 4326, takže je rovnaký ako moje polygóny s PSČ.
point_geo <- st_as_sf (geokódované_adresy,coords = c (x = "lon", y = "lat"),
crs = 4326)
Geopriestorové spojenia so sv
Teraz, keď som nastavil svoje dve množiny údajov, je výpočet sip pre každú adresu pomocou sf jednoduchý st_join () funkcia. Syntax:
st_join (point_sf_object, polygon_sf_object, join = join_type)
V tomto príklade chcem bežať st_join () na geokódovaných bodoch ako prvých a PSČ na polygónoch ako druhých. Je to takzvaný formát ľavého spojenia: Všetky sú zahrnuté body v prvých údajoch (geokódované adresy), ale zodpovedajú iba body v druhom údaji (PSČ). Nakoniec je môj typ spojenia st_within, pretože chcem, aby bol zápas iba body.
my_results <- st_join (point_geo, zipcode_geo,join = st_within)
To je všetko! Teraz, keď sa pozriem na svoje výsledky vytlačením niekoľkých najdôležitejších stĺpcov, uvidíte, že každá adresa má PSČ (v stĺpci „meno“).
print (my_results [, c ("dopyt", "meno", "geometria")])# Jednoduchá zbierka funkcií s 2 funkciami a 2 poliami # typ geometrie: BOD # rozmer: XY # bbox: xmin: -71,39105 ymin: 42,31348 xmax: -71,03673 ymax: 42,34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + údaj = WGS84 + no_defs # geometria názvu dotazu # 1 492 Old Connecticut Path, Framingham, MA 01701 BOD (-71,39105 42,31348) # 2 250 Northern Ave., Boston, MA 02210 BOD (-71,03673 42,34806)
Mapovanie bodov a mnohouholníkov pomocou tmap
Ak chcete mapovať body a mnohouholníky, existuje jeden spôsob, ako to urobiť pomocou tmap:
tm_shape (zipcode_geo) +tm_fill () +
tm_shape (my_results) +
tm_bubbles (col = "red", veľkosť = 0,25)
Snímka obrazovky Sharon Machlis, Chcete viac R tipov? Prejdite na stránku „Robte viac s R“!
