SF를 사용하여 R에서 공간 분석을 수행하는 방법

어디에서 투표합니까? 의원은 누구입니까? 귀하의 우편 번호는 무엇입니까? 이러한 질문에는 지리 공간적으로 공통점이 있습니다. 대답은 점이 속하는 다각형을 결정하는 것입니다.

이러한 계산은 종종 특수 GIS 소프트웨어로 수행됩니다. 하지만 R에서도 쉽게 할 수 있습니다. 세 가지가 필요합니다.

  1. 위도와 경도를 찾기 위해 주소를 지오 코딩하는 방법 
  2. 우편 번호 다각형 경계를 설명하는 셰이프 파일 과 
  3. sf 패키지.

지오 코딩의 경우 일반적으로 geocod.io API를 사용합니다. 하루에 2,500 회 조회가 무료이며 멋진 R 패키지가 있지만 사용하려면 (무료) API 키가 필요합니다. 이 기사의 복잡성을 해결하기 위해 무료 오픈 소스 Open Street Map Nominatim API를 사용할 것입니다. 키가 필요하지 않습니다. 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., Boston, MA")

geocode_OSM을 사용하면 지오 코딩이 간단합니다. 위도 및 경도를 포함하여 처음 세 열을 인쇄하여 결과를 볼 수 있습니다.

geocoded_addresses <-geocode_OSM (주소)

print (geocoded_addresses [, 1 : 3])

위도 조회

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

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

ZIP 코드 셰이프 파일을 얻는 방법에는 여러 가지가 있습니다. 가장 쉬운 방법은 미국 인구 조사국의 우편 번호 표 영역 일 것입니다. 이는 미국 우편 서비스 경계와 정확히 동일하지는 않지만 유사합니다.

미국 인구 조사국에서 직접 ZCTA 파일을 다운로드 할 수 있지만 국가 전체를위한 파일입니다. 대용량 데이터 파일이 마음에 들지 않는 경우에만 그렇게하십시오. 

단일 주에 대한 ZCTA 파일을 다운로드 할 수있는 한 곳은 Census Reporter입니다. 인구와 같은 주별로 데이터를 검색 한 다음 지역에 우편 번호를 추가하고 데이터를 shapefile로 다운로드를 선택합니다.

다운로드 한 파일의 압축을 수동으로 풀 수는 있지만 R에서는 더 쉽습니다. 여기서는 unzip()다운로드 한 파일에 대해 base R의 기능을 사용하고 ma_zip_shapefile이라는 프로젝트 하위 디렉토리에 압축을 풉니 다. 이 junkpaths = TRUE인수는 zip 파일의 이름을 기반으로 다른 하위 디렉토리를 추가하여 압축을 풀고 싶지 않다고 말합니다.

unzip ( "data / acs2017_5yr_B01003_86000US02648.zip",

exdir = "ma_zip_shapefile", junkpaths = TRUE,

덮어 쓰기 = TRUE)

SF를 사용한 지리 공간 가져 오기 및 분석

이제 마지막으로 지리 공간 작업입니다. sf의 st_read()기능을 사용하여 shapefile을 R로 가져올 것 입니다.

zipcode_geo <-st_read ( "ma_zip_shapefile / acs2017_5yr_B01003_86000US02648.shp") # 데이터 소스`/Users/smachlis/Documents/MoreWithR/ma_zip_shape.shp '를 사용하여`acs2017_5yr_B01003_86000US02648'계층 읽기 기능 및 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는 기본적으로 3 차원 지구 (지구)의 영역을 2 차원 좌표 (위도와 경도)로 변환하는 데 사용 된 시스템을 나타냅니다  . 이것은 많은 다른 좌표 참조 시스템 이 있기 때문에 중요 합니다. 내 우편 번호 다각형과 주소 지점이 동일한 것을 사용하여 제대로 정렬되도록하고 싶습니다.

참고 :이 파일에는 매사추세츠 주 전체에 대한 다각형이 포함되어 있지만 필요하지 않습니다. 따라서 매사추세츠 행을 필터링하여

zipcode_geo <-dplyr :: filter (zipcode_geo,

name! = "Massachusetts")

tmap을 사용하여 shapefile 매핑

다각형 데이터를 매핑 할 필요는 없지만 형상 파일이 내가 기대하는 것과 같은지 확인하는 것이 좋습니다. tmap의 qtm()(빠른 테마 맵의 약자) 기능 을 사용하여 sf 객체의 빠른 플롯을 수행 할 수 있습니다 .

qtm (zipcode_geo) +

tm_legend (show = FALSE)

Sharon Machlis가 촬영 한 스크린,

그리고 실제로 우편 번호가 될 수있는 다각형이있는 매사추세츠 기하학이있는 것 같습니다.

다음으로 지오 코딩 된 주소 데이터를 사용하고 싶습니다. 이것은 현재 일반 데이터 프레임이지만 올바른 좌표계를 사용하는 SF 지리 공간 객체로 변환해야합니다.

sf의 st_as_sf()기능으로 할 수 있습니다 . (참고 : 공간 데이터에서 작동하는 sf 패키지 함수는 st_"공간"및 "시간적"을 나타내는로 시작 합니다.)

st_as_sf()몇 가지 인수가 필요합니다. 아래 코드에서 첫 번째 인수는 변환 할 객체 인 내 지오 코딩 된 주소입니다. 두 번째 인수 벡터는 x (경도) 및 y (위도) 값이있는 열을 함수에 알려줍니다. 세 번째는 좌표 참조 시스템을 4326으로 ​​설정하므로 내 우편 번호 다각형과 동일합니다.

point_geo <-st_as_sf (geocoded_addresses,

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

crs = 4326)

SF와 지리 공간 결합

이제 두 개의 데이터 세트를 설정 했으므로 sf의 st_join()기능을 사용하여 각 주소의 우편 번호를 쉽게 계산할 수 있습니다. 구문 :

st_join (point_sf_object, polygon_sf_object, join = join_type)

이 예에서는 st_join()먼저 지오 코딩 된 지점 에서 실행 하고 두 번째로 우편 번호 다각형 에서 실행하려고합니다 . 소위 왼쪽 조인 형식입니다. 첫 번째 데이터 (지오 코딩 된 주소)의 모든 포인트가 포함되지만 일치하는 두 번째 (우편 번호) 데이터의 포인트 만 포함됩니다. 마지막으로, 내 조인 유형은 st_within입니다. 매치가 포인트가되기를 원하기 때문입니다. 

my_results <-st_join (point_geo, zipcode_geo,

조인 = st_within)

그게 다야! 이제 가장 중요한 열 몇 개를 인쇄하여 결과를 보면 각 주소에 우편 번호가있는 것을 볼 수 있습니다 ( "이름"열에 있음). 

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 # 쿼리 이름 geometry # 1492 Old Connecticut Path, Framingham, MA 01701 POINT (-71.39105 42.31348) # 2250 Northern Ave., Boston, MA 02210 POINT (-71.03673 42.34806)

tmap으로 점과 다각형 매핑

포인트와 폴리곤을 매핑하려면 tmap을 사용하는 한 가지 방법이 있습니다.

tm_shape (zipcode_geo) +

tm_fill () +

tm_shape (my_results) +

tm_bubbles (col = "red", 크기 = 0.25)

Sharon Machlis의 스크린 샷,

더 많은 R 팁을 원하십니까? "R로 더 많은 작업 수행"페이지로 이동하십시오!