Jak przeprowadzić analizę przestrzenną w R z sf

Gdzie głosujesz? Kim jesteście ustawodawcami? Jaki jest twój kod pocztowy? Te pytania mają coś wspólnego pod względem geoprzestrzennym: odpowiedź polega na określeniu, do którego wielokąta należy dany punkt.

Takie obliczenia są często wykonywane za pomocą specjalistycznego oprogramowania GIS. Ale są też łatwe do zrobienia w R. Potrzebujesz trzech rzeczy:

  1. Sposób geokodowania adresów w celu znalezienia szerokości i długości geograficznej; 
  2. Pliki kształtów określające granice wielokątów kodu pocztowego; i 
  3. Pakiet sf.

Do geokodowania zazwyczaj używam API geocod.io. Jest darmowy dla 2500 wyszukiwań dziennie i ma ładny pakiet R, ale potrzebujesz (darmowego) klucza API, aby z niego korzystać. Aby obejść tę złożoność tego artykułu, skorzystam z bezpłatnego, otwartego źródła API Open Street Map Nominatim. Nie wymaga klucza. Pakiet tmaptools ma funkcję, geocode_OSM()używającą tego API.

Importowanie i przygotowywanie danych geoprzestrzennych

Będę używać pakietów sf, tmaptools, tmap i dplyr. Jeśli chcesz kontynuować, załaduj każdy z nich pacman::p_load()lub zainstaluj go, którego jeszcze nie ma w systemie install.packages(), a następnie załaduj każdy z library().

W tym przykładzie utworzę wektor z dwoma adresami, naszym biurem w Framingham w stanie Massachusetts i biurem RStudio w Bostonie.

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

„250 Northern Ave., Boston, MA”)

Geokodowanie jest proste dzięki geocode_OSM. Możesz zobaczyć wyniki, drukując pierwsze trzy kolumny, w tym szerokość i długość geograficzną:

geocoded_addresses <- geocode_OSM (adresy)

drukuj (geokodowane_adresy [, 1: 3])

zapytanie lat lon

# 1492 Old Connecticut Path, Framingham, MA 42.31348 -71.39105

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

Istnieje kilka sposobów na pobranie plików kształtu kodu pocztowego. Najłatwiejsze są prawdopodobnie obszary tabel kodów pocztowych US Census Bureau, które są podobne, jeśli nie dokładnie takie same, jak granice usług pocztowych USA.

Możesz pobrać plik ZCTA bezpośrednio z US Census Bureau, ale jest to plik dla całego kraju. Zrób to tylko wtedy, gdy nie przeszkadza Ci duży plik danych. 

Jednym miejscem do pobrania pliku ZCTA dla pojedynczego stanu jest Census Reporter. Wyszukaj dowolne dane według stanu, takie jak populacja, a następnie dodaj kod pocztowy do obszaru geograficznego i wybierz pobieranie danych jako plik kształtu.

Mógłbym ręcznie rozpakować pobrany plik, ale jest to łatwiejsze w R. Tutaj używam podstawowej unzip()funkcji R na pobranym pliku i rozpakowuję go do podkatalogu projektu o nazwie ma_zip_shapefile. Ten junkpaths = TRUEargument mówi, że nie chcę rozpakowywać, dodając kolejny podkatalog na podstawie nazwy pliku zip.

unzip ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

overwrite = TRUE)

Import i analiza danych geoprzestrzennych za pomocą sf

Teraz wreszcie trochę pracy geoprzestrzennej. Importuję shapefile do R używając st_read()funkcji sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Odczytywanie warstwy `acs2017_5yr_B01003_86000US02648 'ze źródła danych` /Users/smachlis/Documents/MoreWip.s_000_008z' 'ES486Prosty_BlZ_006' cechy i 4 pola # typ geometrii: MULTIPOLYGON # wymiar: XY # bbox: xmin: -73.50821 ymin: 41.18705 xmax: -69.85886 ymax: 42.95774 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs

st_read()Dołączam odpowiedź konsoli podczas uruchamiania, ponieważ są tam wyświetlane informacje: plik epsg. To mówi, jaki system odniesienia za pomocą współrzędnych został użyty do utworzenia pliku . Tutaj było 4326. Bez wchodzenia zbyt głęboko w chwasty, epsg zasadniczo wskazuje,  jaki system został użyty do przetłumaczenia obszarów na trójwymiarowej kuli ziemskiej - Ziemi - na współrzędne dwuwymiarowe (szerokość i długość) . Jest to ważne, ponieważ istnieje wiele różnych układów odniesienia za pomocą współrzędnych. Chcę, aby moje wielokąty kodu pocztowego i punkty adresowe używały tego samego, aby były prawidłowo wyrównane.

Uwaga: zdarza się, że ten plik zawiera wielokąt dla całego stanu Massachusetts, którego nie potrzebuję. Więc odfiltruję ten wiersz Massachusetts z

zipcode_geo <- dplyr :: filter (zipcode_geo,

name! = „Massachusetts”)

Mapowanie shapefile za pomocą tmap

Mapowanie danych wielokątów nie jest konieczne, ale dobrze jest sprawdzić mój plik shapefile, aby zobaczyć, czy geometria jest tym, czego oczekuję. Możesz zrobić szybki wykres obiektu sf za pomocą funkcji tmap qtm()(skrót od quick theme map).

qtm (zipcode_geo) +

tm_legend (show = FALSE)

Ekrany nakręcone przez Sharon Machlis,

I wygląda na to, że rzeczywiście mam geometrię Massachusetts z wielokątami, które mogą być kodami pocztowymi.

Następnie chcę użyć geokodowanych danych adresowych. Obecnie jest to zwykła ramka danych, ale musi zostać przekonwertowana na obiekt geoprzestrzenny sf z odpowiednim układem współrzędnych.

Możemy to zrobić za pomocą st_as_sf()funkcji sf . (Uwaga: funkcje pakietu sf, które działają na danych przestrzennych, zaczynają się od st_, co oznacza „przestrzenny” i „czasowy”).

st_as_sf()przyjmuje kilka argumentów. W poniższym kodzie pierwszym argumentem jest obiekt do przekształcenia - moje geokodowane adresy. Wektor drugiego argumentu informuje funkcję, które kolumny mają wartości x (długość) i y (szerokość geograficzna). Trzeci ustawia system odniesienia za pomocą współrzędnych na 4326, więc jest taki sam jak wielokąty mojego kodu pocztowego.

point_geo <- st_as_sf (geocoded_addresses,

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

crs = 4326)

Połączenia geoprzestrzenne z sf

Teraz, gdy skonfigurowałem moje dwa zestawy danych, obliczenie kodu pocztowego dla każdego adresu jest łatwe dzięki st_join()funkcji sf . Składnia:

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")])

# Prosty zbiór funkcji z 2 funkcjami i 2 polami # typ geometrii: POINT # wymiar: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # nazwa zapytania geometria # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

Mapowanie punktów i wielokątów za pomocą tmap

Jeśli chcesz zmapować punkty i wielokąty, oto jeden ze sposobów na zrobienie tego za pomocą tmap:

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (my_results) +

tm_bubbles (col = "red", size = 0,25)

Zrzut ekranu: Sharon Machlis,

Chcesz więcej wskazówek R? Przejdź na stronę „Zrób więcej z R”!