Как да направя пространствен анализ в R с sf

Къде гласувате? Кои сте законодателите? Какъв е вашият пощенски код? Тези въпроси имат нещо общо геоспатично: Отговорът включва определяне в кой полигон попада една точка.

Такива изчисления често се правят със специализиран ГИС софтуер. Но те също са лесни за изпълнение в R. Имате нужда от три неща:

  1. Начин за геокодиране на адреси за намиране на географска ширина и дължина; 
  2. Shapefiles, които очертават границите на полигона на пощенския код; и 
  3. Пакетът sf.

За геокодиране обикновено използвам API на geocod.io. Безплатно е за 2500 търсения на ден и има хубав R пакет, но за да го използвате, ви е необходим (безплатен) API ключ. За да заобиколя тази част от сложността на тази статия, ще използвам безплатния API с отворени кодове Open Street Nominatim. Не изисква ключ. Пакетът tmaptools има функция geocode_OSM(), за да използва този API.

Импортиране и подготовка на геопространствени данни

Ще използвам пакетите sf, tmaptools, tmap и dplyr. Ако искате да продължите, заредете всеки с pacman::p_load()или инсталирайте още, който все още не е във вашата система install.packages(), след това заредете всеки с library().

За този пример ще създам вектор с два адреса, офиса ни във Фрамингъм, Масачузетс и офиса RStudio в Бостън.

адреси <- c ("492 Old Connecticut Path, Framingham, MA",

"250 Northern Ave., Бостън, Масачузетс")

Геокодирането е лесно с geocode_OSM. Можете да видите резултатите, като отпечатате първите три колони, включително географска ширина и дължина:

geocoded_addresses <- geocode_OSM (адреси)

печат (геокодирани_адреси [, 1: 3])

заявка lat lon

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

# 2 250 Northern Ave., Бостън, Масачузетс 42.34806 -71.03673

Има няколко начина за получаване на шейпфайлове с пощенски код. Най-лесният е може би зоните за табулация на пощенския код на Бюрото за преброяване на населението на САЩ, които са подобни, ако не и съвсем същите като границите на пощенската служба на САЩ.

Можете да изтеглите ZCTA файл директно от Бюрото за преброяване на населението на САЩ, но това е файл за цялата страна. Правете това само ако нямате нищо против голям файл с данни. 

Едно място за изтегляне на ZCTA файл за една държава е Census Reporter. Потърсете всякакви данни по щати, като население, и след това добавете пощенски код към географията и изберете изтегляне на данни като шейпфайл.

Бих могъл да разархивирам изтегления си файл ръчно, но в R. е по-лесно. Тук използвам unzip()функцията base R на изтеглен файл и го разархивирам в поддиректория на проекта с име ma_zip_shapefile. Този junkpaths = TRUEаргумент казва, че не искам да разархивирам добавяне на друга поддиректория въз основа на името на zip файла.

разархивирайте ("data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

презаписване = TRUE)

Геопространствен внос и анализ със sf

Сега накрая някаква геопространствена работа. Ще импортирам шейп файла в R, използвайки st_read()функцията на sf .

zipcode_geo <- st_read ("ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # Четене на слоя `acs2017_5yr_B01003_86000US02648 'от източник на данни` /Users/smachlis/Documents/MoreWithRsha0_01 характеристики и 4 полета # тип геометрия: MULTIPOLYGON # измерение: 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()Включих отговора на конзолата при стартиране, защото там се показва някаква информация: epsg. Това казва каква координатна система за справка е използвана за създаване на файла . Тук беше 4326. Без да навлиза прекалено дълбоко в плевелите, epsg основно показва  коя система е била използвана за превод на области на триизмерен глобус - Земята - в двуизмерни координати (географска ширина и дължина) . Това е важно, защото има много различни координатни референтни системи. Искам моите полигони с пощенски код и адресните точки да използват един и същ, така че да се подреждат правилно.

Забележка: Случайно този файл включва полигон за целия щат Масачузетс, който не ми трябва. Така че ще филтрирам този ред в Масачузетс

zipcode_geo <- dplyr :: filter (zipcode_geo,

name! = "Масачузетс")

Съставяне на шейп файла с tmap

Картирането на данните на полигона не е необходимо, но е хубава проверка на моя шейпфайл, за да се види дали геометрията е това, което очаквам. Можете да направите бърз парцел на sf обект с функцията на tmap qtm()(съкратено от бърза карта на теми).

qtm (zipcode_geo) +

tm_legend (show = FALSE)

Екрани, заснети от Шарън Маклис,

И изглежда, че наистина имам геометрия на Масачузетс с полигони, които могат да бъдат пощенски кодове.

След това искам да използвам геокодирани данни за адрес. Понастоящем това е обикновен кадър от данни, но той трябва да бъде преобразуван в sf геопространствен обект с правилната координатна система.

Можем да направим това с st_as_sf()функцията на sf . (Забележка: sf функциите на пакета, които работят върху пространствени данни, започват с st_, което означава „пространствен“ и „времеви“.)

st_as_sf()взема няколко аргумента. В кода по-долу първият аргумент е обектът за трансформиране - моите геокодирани адреси. Вторият вектор на аргумента казва на функцията кои колони имат стойностите x (дължина) и y (географска ширина). Третият задава координатната референтна система на 4326, така че тя е същата като моите полигони с пощенски код.

point_geo <- st_as_sf (геокодирани_адреси,

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

crs = 4326)

Геопространствени съединения със sf

Сега, когато настроих двата си набора от данни, изчисляването на пощенския код за всеки адрес е лесно с st_join()функцията на sf . Синтаксисът:

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

# Проста колекция от функции с 2 функции и 2 полета # тип геометрия: POINT # измерение: XY # bbox: xmin: -71.39105 ymin: 42.31348 xmax: -71.03673 ymax: 42.34806 # epsg (SRID): 4326 # proj4string: + proj = longlat + datum = WGS84 + no_defs # геометрия на името на заявката # 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)

Картографски точки и полигони с tmap

Ако искате да картографирате точките и полигоните, ето един начин да го направите с tmap:

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (моите_резултати) +

tm_bubbles (col = "червен", размер = 0,25)

Снимка на екрана от Шарън Маклис,

Искате още R съвети? Насочете се към страницата „Направете повече с R“!