Kako napraviti prostornu analizu u R sa sf

Gdje glasate? Tko ste vi zakonodavci? Koji je vaš poštanski broj? Ova pitanja imaju nešto zajedničko u geospatiji: Odgovor uključuje određivanje u koji poligon spada točka.

Takvi se izračuni često rade sa specijaliziranim GIS softverom. Ali to je lako učiniti i u R. Potrebne su vam tri stvari:

  1. Način geokodiranja adresa za pronalaženje zemljopisne širine i dužine; 
  2. Datoteke oblika koje ocrtavaju granice poligona poštanskog broja; i 
  3. SF paket.

Za geokodiranje obično koristim API geocod.io. Besplatno je za 2.500 pregleda dnevno i ima lijep R paket, ali za njegovo korištenje potreban vam je (besplatni) API ključ. Da bih zaobišao malo složenosti ovog članka, poslužit ću se besplatnim API-jem otvorenih mapa Open Map Nominatim. Ne treba ključ. Paket tmaptools ima funkciju,, geocode_OSM()da koristi taj API.

Uvoz i priprema geoprostornih podataka

Koristit ću pakete sf, tmaptools, tmap i dplyr. Ako želite nastaviti, učitajte svaki sa pacman::p_load()ili instalirajte bilo koji koji još nije na vašem sustavu install.packages(), a zatim učitajte svaki sa library().

Za ovaj primjer stvorit ću vektor s dvije adrese, našim uredom u Framinghamu u Massachusettsu i uredom RStudio u Bostonu.

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

"250 Northern Ave., Boston, MA")

Geokodiranje je jednostavno s geocode_OSM. Rezultate možete vidjeti ispisom prva tri stupca, uključujući zemljopisnu širinu i dužinu:

geokodirane_adrese <- geocode_OSM (adrese)

ispis (geocoded_addresses [, 1: 3])

upit lat lon

# 1 492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105

# 2 250 Northern Ave., Boston, MA 42.34806 -71.03673

Postoji nekoliko načina za dobivanje datoteka u obliku poštanskog broja. Najlakše je vjerojatno područje tabuliranja poštanskog broja američkog Zavoda za popis stanovništva, koje je slično ako ne i potpuno isto kao granice američkih poštanskih službi.

ZCTA datoteku možete preuzeti izravno s američkog ureda za popis stanovništva, ali to je datoteka za cijelu zemlju. Učinite to samo ako vam ne smeta velika podatkovna datoteka. 

Census Reporter je jedno mjesto za preuzimanje ZCTA datoteke za jednu državu. Potražite bilo koje podatke prema državi, poput stanovništva, a zatim dodajte poštanski broj u zemljopis i odaberite podatke za preuzimanje kao datoteku oblika.

Mogao bih ručno raspakirati preuzetu datoteku, ali u R. je lakše. Ovdje koristim osnovnu unzip()funkciju R na preuzetoj datoteci i raspakiram je u poddirektorij projekta pod nazivom ma_zip_shapefile. Taj junkpaths = TRUEargument kaže da ne želim raspakirati dodajući još jedan poddirektorij na temelju imena zip datoteke.

raspakirajte ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

prepiši = TRUE)

Geoprostorni uvoz i analiza sf

Sad napokon neki geoprostorni rad. Datoteku oblika uvest ću u R pomoću st_read()funkcije sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Sloj za čitanje `acs2017_5yr_B01003_86000US02648 'iz izvora podataka` /Users/smachlis/Documents/MoreWithR01_01' značajke i 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_defs

Prilikom izvođenja uključio sam odgovor konzole st_read()jer su tamo prikazane neke informacije: epsg. To govori koji je koordinatni referentni sustav korišten za stvaranje datoteke . Ovdje je bilo 4326. Ne ulazeći preduboko u korov, epsg u osnovi pokazuje  koji je sustav korišten za prevođenje područja na trodimenzionalnoj globusu - Zemlji - u dvodimenzionalne koordinate (zemljopisnu širinu i dužinu) . To je važno jer postoji puno različitih koordinatnih referentnih sustava. Želim da moji poligoni poštanskih brojeva i adresne točke koriste isti, tako da se pravilno poredaju.

Napomena: Ova datoteka slučajno sadrži poligon za cijelu državu Massachusetts, koji mi nije potreban. Pa ću filtrirati taj Massachusettsov red s

zipcode_geo <- dplyr :: filter (zipcode_geo,

ime! = "Massachusetts")

Mapiranje shapefile-a s tmap

Mapiranje podataka poligona nije potrebno, ali lijepa je provjera moje shapefile datoteke da vidim je li geometrija ono što očekujem. Možete brzo zacrtati sf objekt pomoću funkcije tmap qtm()(kratica za brzu mapu tema).

qtm (zipcode_geo) +

tm_legend (show = FALSE)

Ekrane snimila Sharon Machlis,

Izgleda da zaista imam geometriju Massachusettsa s poligonima koji bi mogli biti poštanski brojevi.

Dalje želim koristiti podatke geokodirane adrese. Ovo je trenutno običan podatkovni okvir, ali ga treba pretvoriti u sf geoprostorni objekt s pravim koordinatnim sustavom.

To možemo učiniti sf-ovom st_as_sf()funkcijom. (Napomena: sf funkcije paketa koje djeluju na prostornim podacima započinju s st_, što je skraćenica za "prostorni" i "vremenski".)

st_as_sf()uzima nekoliko argumenata. U donjem kodu prvi je argument objekt koji treba transformirati - moje geokodirane adrese. Drugi vektor argumenata govori funkciji koji stupci imaju vrijednosti x (zemljopisna dužina) i y (zemljopisna širina). Treća postavlja referentni koordinatni sustav na 4326, pa je to isto kao i moji poligoni za poštanski broj.

point_geo <- st_as_sf (geokodirane_adrese,

koordinate = c (x = "dugo", y = "lat"),

crs = 4326)

Geoprostorni spojevi sf

Sad kad sam postavio svoja dva skupa podataka, izračunavanje poštanskog broja za svaku adresu lako je pomoću st_join()funkcije sf . Sintaksa:

st_join (point_sf_object, polygon_sf_object, join = join_type)

In this example, I want to run st_join() on the geocoded points first and the ZIP code polygons second. It’s a so-called left join format: All points in the first data (geocoded addresses) are included, but only points in the second (ZIP code) data that match. Finally, my join type is st_within, since I want the match to be points within. 

my_results <- st_join(point_geo, zipcode_geo,

join = st_within)

That’s it! Now if I look at my results by printing out several of the most important columns, you”ll see each address has a ZIP code (in the “name” column). 

print(my_results[,c("query", "name", "geometry")])

# Jednostavna kolekcija značajki s 2 značajke i 2 polja # vrsta geometrije: TOČKA # dimenzija: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # geometrija imena upita # 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)

Mapiranje točaka i poligona s tmapom

Ako želite mapirati točke i poligone, evo jednog načina da to učinite pomoću tmap:

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (moji_rezultati) +

tm_bubbles (col = "crvena", veličina = 0,25)

Snimak ekrana Sharon Machlis,

Želite još R savjeta? Krenite na stranicu "Učini više s R"!