Kje voliš? Kdo ste zakonodajalci? Kakšna je vaša poštna številka? Ta vprašanja imajo nekaj skupnega z geospatijo: odgovor vključuje določitev, v kateri poligon spada točka.
Takšni izračuni se pogosto izvajajo s specializirano GIS programsko opremo. Toda v R. jih je enostavno narediti. Potrebujete tri stvari:
- Način za geokodiranje naslovov za iskanje zemljepisne širine in dolžine;
- Datoteke v obliki, ki orisajo meje poligona poštne številke; in
- Sf paket.
Za geokodiranje običajno uporabljam API geocod.io. Brezplačno je 2.500 iskanj na dan in ima lep paket R, vendar za njegovo uporabo potrebujete (brezplačen) ključ API. Da bi se izognil tej zapletenosti tega članka, bom uporabil brezplačen odprtokodni API Open Map Nominatim. Ne zahteva ključa. Paket tmaptools ima funkcijo, geocode_OSM (), da uporabite ta API.
Uvoz in priprava geoprostorskih podatkov
Uporabljal bom pakete sf, tmaptools, tmap in dplyr. Če želite nadaljevati, naložite vsakega z pacman :: p_load () ali namestite katerega koli, ki še ni v vašem sistemu z install.packages (), nato naložite vsakega z knjižnica().
Za ta primer bom ustvaril vektor z dvema naslovoma, našo pisarno v Framinghamu v Massachusettsu in pisarno RStudio v Bostonu.
naslovi <- c ("492 Old Connecticut Path, Framingham, MA","250 Northern Ave., Boston, MA")
Geocode je enostavno z geocode_OSM. Rezultate si lahko ogledate tako, da natisnete prve tri stolpce, vključno z zemljepisno širino in dolžino:
geocoded_addresses <- geocode_OSM (naslovi)tisk (geododirani_naslovi [, 1: 3])
poizvedba lat lon
# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105
# 2 250 Northern Ave., Boston, MA 42.34806 -71.03673
Obstaja več načinov za pridobivanje datotek oblike poštne številke. Najlažje je verjetno območje za tabelariranje poštnih številk urada za popis prebivalstva ZDA, ki je podobno, če ne povsem enako mejam ameriške poštne službe.
Datoteko ZCTA lahko prenesete neposredno iz urada za popis prebivalstva ZDA, vendar je datoteka za celotno državo. To storite le, če vas ne moti velika podatkovna datoteka.
Eno mesto za prenos datoteke ZCTA za posamezno državo je Census Reporter. Poiščite podatke po zvezni državi, na primer prebivalstvo, nato geografski lokaciji dodajte poštno številko in izberite datoteko za prenos kot datoteko oblike.
Svojo preneseno datoteko bi lahko razpakiral ročno, vendar je v R. lažje. Tu uporabljam osnovno R razpakiraj () funkcijo v preneseni datoteki in jo razpakirajte v podimenik projekta z imenom ma_zip_shapefile. To junkpaths = TRUE argument pravi, da ne želim razpakirati dodajanja novega podimenika na podlagi imena zip datoteke.
razpakirajte ("data / acs2017_5yr_B01003_86000US02648.zip",exdir = "ma_zip_shapefile", junkpaths = TRUE,
prepiši = TRUE)
Geoprostorski uvoz in analiza s sf
Zdaj končno nekaj geoprostorskega dela. Datoteko oblike bom uvozil v R s pomočjo sf-jev st_read () funkcijo.
zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Bralni sloj `acs2017_5yr_B01003_86000US02648 'iz vira podatkov` /Users/smachlis/Documents/MoreWithR01_01' funkcije in 4 polja # vrsta geometrije: MULTIPOLYGON # dimenzija: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defsMed izvajanjem sem vključil odgovor konzole st_read () ker je tam prikazanih nekaj informacij: epsg. To pravi kakšen koordinatni referenčni sistem je bil uporabljen za izdelavo datoteke. Tu je bilo 4326. Epsg v bistvu kaže, ne da bi se preveč zapletal v plevelkateri sistem je bil uporabljen za prevajanje območij na tridimenzionalnem globusu - Zemlji - v dvodimenzionalne koordinate (zemljepisne širine in dolžine). To je pomembno, ker obstajajo veliko različnih koordinatnih referenčnih sistemov. Želim, da bi moji poligoni poštne številke in naslovne točke uporabljali istega, zato se pravilno ujemajo.
Opomba: Ta datoteka vsebuje poligon za celotno zvezno državo Massachusetts, ki ga ne rabim. Torej bom filtriral tisto vrstico v Massachusettsu
zipcode_geo <- dplyr :: filter (zipcode_geo,name! = "Massachusetts")
Preslikava datoteke oblike z datoteko tmap
Kartiranje podatkov poligona ni potrebno, je pa lep pregled moje datoteke oblike, da vidim, ali je geometrija takšna, kot pričakujem. Hitro lahko narišete sf objekt s tmapi qtm () (kratica za hitri zemljevid tem).
qtm (zipcode_geo) +tm_legend (show = FALSE)
Zaslone je posnela Sharon Machlis, In zdi se, da resnično imam geometrijo Massachusettsa s poligoni, ki bi lahko bile poštne številke.
Nato želim uporabiti podatke o geokodiranem naslovu. To je trenutno navaden podatkovni okvir, vendar ga je treba pretvoriti v sf geoprostorski objekt s pravim koordinatnim sistemom.
To lahko storimo s sf-ji st_as_sf () funkcijo. (Opomba: funkcije paketa sf, ki delujejo na prostorskih podatkih, se začnejo z st_, kar pomeni "prostorski" in "časovni.")
st_as_sf () vzame več argumentov. V spodnji kodi je prvi argument predmet za pretvorbo - moji geokodirani naslovi. Drugi argument argumenta pove funkciji, kateri stolpci imajo vrednosti x (zemljepisna dolžina) in y (zemljepisna širina). Tretji nastavi koordinatni referenčni sistem na 4326, torej je enak kot moji poligoni za poštno številko.
point_geo <- st_as_sf (geokodirani_naslovi,coords = c (x = "lon", y = "lat"),
crs = 4326)
Geoprostorski spoji s sf
Zdaj, ko sem nastavil svoja dva nabora podatkov, je sf-jevo izračunavanje poštne številke za vsak naslov enostavno st_join () funkcijo. Sintaksa:
st_join (point_sf_object, polygon_sf_object, join = join_type)
V tem primeru želim teči st_join () najprej na geokodirane točke, drugi pa na poligone poštne številke. Gre za tako imenovano obliko levega združevanja: Vse vključene so točke v prvih podatkih (geokodirani naslovi), vendar le točke v podatkih drugega (poštna številka), ki se ujemajo. Končno je moj tip pridružitve st_within, ker želim, da je tekma znotraj točk.
moji_rezultati <- st_join (point_geo, zipcode_geo,pridruži = st_within)
To je to! Če pogledam svoje rezultate tako, da natisnem več najpomembnejših stolpcev, boste videli, da ima vsak naslov poštno številko (v stolpcu »ime«).
natisni (moji_rezultati [, c ("poizvedba", "ime", "geometrija")))# Enostavna zbirka funkcij z dvema značilnostma in dvema poljem + datum = WGS84 + no_defs # geometrija imena poizvedbe # 1 492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2 250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)
Preslikava točk in poligonov s tmap
Če želite zemljevide točk in mnogokotnikov preslikati na naslednji način, lahko to naredite s tmap:
tm_shape (zipcode_geo) +tm_fill () +
tm_shape (moji_rezultati) +
tm_bubbles (col = "rdeča", velikost = 0,25)
Posnetek zaslona Sharon Machlis, Želite več nasvetov R? Pojdite na stran »Naredi več z R«!
