Ohjelmointi

Kuinka tehdä alueellinen analyysi R: ssä sf: n kanssa

Missä äänestät? Keitä te olette lainsäätäjiä? Mikä on postinumero? Näillä kysymyksillä on jotain geospatiaalisesti yhteistä: Vastaus sisältää sen määrittämisen, mihin polygoniin piste kuuluu.

Tällaiset laskelmat tehdään usein erikoistuneilla paikkatietojärjestelmillä. Mutta ne on myös helppo tehdä R.: ssä. Tarvitset kolme asiaa:

  1. Tapa geokoodata osoitteita leveys- ja pituuspiirien löytämiseksi;
  2. Muototiedostot, jotka hahmottavat postinumeron monikulmion rajat; ja
  3. SF-paketti.

Geokoodauksessa käytän yleensä geocod.io-sovellusliittymää. Se on ilmainen 2500 hakua päivässä ja siinä on mukava R-paketti, mutta tarvitset (ilmaisen) API-avaimen sen käyttämiseen. Tämän artikkelin monimutkaisuuden kiertämiseksi käytän ilmaista, avoimen lähdekoodin avoimen katukartan Nominatim-sovellusliittymää. Se ei vaadi avainta. Tmaptools-paketilla on toiminto, geokoodi_OSM (), käyttää kyseistä sovellusliittymää.

Paikkatietojen tuominen ja valmistelu

Käytän sf-, tmaptools-, tmap- ja dplyr-paketteja. Jos haluat seurata mukana, lataa kukin pacman :: p_load () tai asenna sellaisia, joita ei vielä ole järjestelmässäsi install.packages ()ja lataa sitten jokaisella kirjasto().

Tässä esimerkissä luon vektorin, jolla on kaksi osoitetta, toimistomme Framinghamissa Massachusettsissa ja RStudion toimisto Bostonissa.

osoitteet <- c ("492 Old Connecticut Path, Framingham, MA",

"250 Northern Ave., Boston, MA")

Geokoodaus on suoraviivaista geokoodin_OSM kanssa. Näet tulokset tulostamalla kolme ensimmäistä saraketta, mukaan lukien leveys- ja pituusaste:

geokoodatut osoitteet <- geokoodin_OSM (osoitteet)

tulosta (geokoodatut_osoitteet [, 1: 3])

kysely lat lon

# 1492 Old Connecticut Path, Framingham, MA 42.31348 -71,39105

# 2250 Northern Ave., Boston, MA 42.34806 -71.03673

Postinumero-muotoisia tiedostoja on useita tapoja saada. Helpoin on luultavasti Yhdysvaltain väestönlaskennatoimiston postinumerotaulukkoalueet, jotka ovat samanlaisia ​​kuin ellei tarkalleen samat kuin Yhdysvaltain postipalvelun rajat.

Voit ladata ZCTA-tiedoston suoraan Yhdysvaltain väestönlaskentatoimistosta, mutta se on koko maata koskeva tiedosto. Tee se vain, jos et välitä suuresta datatiedostosta.

Yksi paikka ladata ZCTA-tiedosto yhdelle osavaltiolle on Census Reporter. Hae mitä tahansa tietoja valtion, kuten väestön, mukaan, lisää sitten postinumero maantieteelliseen alueeseen ja valitse lataustiedot muototiedostoksi.

Voin purkaa ladatun tiedoston manuaalisesti, mutta se on helpompaa R. pura () toiminto ladatulle tiedostolle ja pura se projektin alihakemistoon nimeltä ma_zip_shapefile. Että romupolut = TOSI väite sanoo, etten halua purkaa uuden alihakemiston lisäämistä zip-tiedoston nimen perusteella.

pura ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", roskapolut = TOSI,

korvaa = TOSI)

Paikkatietojen tuonti ja analyysi sf: n kanssa

Nyt vihdoin geospatiaalista työtä. Tuon shapefile-tiedoston R: ään sf: n avulla st_read () toiminto.

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Lukukerros "acs2017_5yr_B01003_86000US02648" tietolähteestä "/Users/smachlis/Documents/MoreWits_ap_Play_app_app_app_app_app_app_app_app_app_ap__set_app_et_ap__set_app_et_painos} ominaisuudet ja 4 kenttää # geometriatyyppi: MULTIPOLYGON # ulottuvuus: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_def

Olen sisällyttänyt konsolin vastauksen käynnissä st_read () koska siellä on joitain tietoja: epsg. Siinä lukee mitä koordinaattien referenssijärjestelmää käytettiin tiedoston luomiseen. Tässä se oli 4326. Epsg osoittaa periaatteessa menemättä liian syvälle rikkaruohoihinmitä järjestelmää käytettiin kolmiulotteisen maapallon - maapallon - alueiden kääntämiseen kaksiulotteisiksi koordinaateiksi (leveys- ja pituusasteet). Tämä on tärkeää, koska on olemassa paljon eri koordinaattijärjestelmiä. Haluan, että postinumeron polygonit ja osoitepisteet käyttävät samaa, joten ne asettuvat riviin oikein.

Huomaa: Tämä tiedosto sisältää satunnaisesti monikulmion koko Massachusettsin osavaltiolle, jota en tarvitse. Joten suodatan pois Massachusettsin rivin

postinumero postinumero_geo <- dplyr :: suodatin (postinumero postinumero,

nimi! = "Massachusetts")

Muotoilutiedoston kartoitus tmapilla

Monikulmion tietojen kartoittaminen ei ole välttämätöntä, mutta se on mukava tarkistaa muotoilutiedostoni nähdäksesi, onko geometria sitä mitä odotan. Voit tehdä pikakuvion sf-objektista tmap: illa qtm () (lyhyt pikakartta).

qtm (postinumero_geo) +

tm_legend (näytä = EPÄTOSI)

Näyttökuvat: Sharon Machlis,

Ja näyttää siltä, ​​että minulla on todellakin Massachusettsin geometria, jossa on polygoneja, jotka voivat olla postinumeroita.

Seuraavaksi haluan käyttää geokoodattuja osoitetietoja. Tämä on tällä hetkellä tavallinen tietokehys, mutta se on muunnettava sf-paikkatietokohteeksi, jolla on oikea koordinaattijärjestelmä.

Voimme tehdä sen sf: llä st_as_sf () toiminto. (Huomaa: Paikkatiedoilla toimivat sf-pakettitoiminnot alkavat st_, joka tarkoittaa "spatiaalista" ja "ajallista".)

st_as_sf () ottaa useita argumentteja. Alla olevassa koodissa ensimmäinen argumentti on muunnettava kohde - geokoodatut osoitteet. Toinen argumenttivektori kertoo funktiolle, mitkä sarakkeet ovat x (pituus) ja y (leveys) arvoja. Kolmas asettaa koordinaattien referenssijärjestelmäksi 4326, joten se on sama kuin postinumeron polygonini.

point_geo <- st_as_sf (geokoodatut_osoitteet,

koordinaatit = c (x = "lon", y = "lat"),

crs = 4326)

Geospatial yhdistyy sf

Nyt kun olen määrittänyt kaksi tietojoukkoani, jokaisen osoitteen postinumero on helppo laskea sf: n avulla st_join () toiminto. Syntaksi:

st_join (point_sf_object, polygon_sf_object, join = join_type)

Tässä esimerkissä haluan juosta st_join () geokoodatuista pisteistä ensin ja postinumero-polygoneista toiseksi. Se on niin kutsuttu vasemman liittymän muoto: Kaikki Ensimmäisten tietojen pisteet (geokoodatut osoitteet) sisältyvät, mutta vain toisessa (postinumero) oleviin tietoihin liittyvät pisteet. Lopuksi, liittymistyyppini on st_ sisällä, koska haluan ottelun olevan pisteitä sisällä.

omatulokset <- st_join (point_geo, postinumero_geo,

liittyä = sis. sisällä)

Se siitä! Nyt kun tarkastelen tuloksia tulostamalla useita tärkeimpiä sarakkeita, näet, että jokaisella osoitteella on postinumero ("nimi" -sarakkeessa).

tulosta (omatulokset [, c ("kysely", "nimi", "geometria")])

# Yksinkertainen ominaisuuksien kokoelma, jossa on 2 ominaisuutta ja 2 kenttää # geometriatyyppi: POINT # ulottuvuus: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + peruspiste = WGS84 + no_defs # kyselyn nimen geometria # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Pisteiden ja polygonien kartoitus tmapilla

Jos haluat kartoittaa pisteet ja polygonit, tämä on yksi tapa tehdä se tmap:

tm_shape (postinumero_geo) +

tm_fill () +

tm_shape (omatulokset) +

tm_bubbles (col = "punainen", koko = 0,25)

Näyttökuva: Sharon Machlis,

Haluatko lisää R-vinkkejä? Suuntaa "Tee enemmän R: llä" -sivulle!