저는 R에서 GIS를 처음 사용하고 lat / long 또는 UTM 좌표를 shapefile에 추가하려고합니다. 여기에서 시카고 (City_Boundaries.shp)의 경계를 다운로드했습니다. http://www.cityofchicago.org/city/en/depts/doit/supp_info/gis_data.html
maptools 및 rgeos 라이브러리를로드했습니다.
library(maptools)
library(rgeos)
library(rgdal)
데이터를 R로 가져오고 영역 16T에 대한 UTM 코드를 추가하려고했습니다.
city<- readShapeSpatial("C:/Users/Luke.Gerdes/Documents/GIS Files/Chicago/City_Boundary.shp")
proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")
그러나 결과 데이터는 나에게 의미가 없습니다. "city"내의 "coords"슬롯을 살펴보면 좌표 값은 예를 들어 X = 1092925 및 Y = 1944820입니다. GPS 좌표를 찾기 위해 외부 도구 ( http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html )를 사용했습니다 . 결과 : lat = 17.511, long = -81.42. 이것은 자메이카와 온두라스 해안 사이에 있습니다. 나에게는 shapefile 좌표가 자신의 우주에 존재하는 것 같습니다. shapefile이라는 이름에서 알 수 있듯이 좌표는 도시의 모양을 정확하게 묘사하지만 이러한 좌표는 지구본에 자동으로 매핑되지 않습니다.
내 궁극적 인 목표는 시카고 내에서 여러 지역 태그가 지정된 이벤트가 발생했는지 확인하는 것입니다. 위도 / 경도로 나열된 포인트를 아래에 따라 UTM으로 변환하는 것이 편안합니다.
SP <- SpatialPoints(cbind(-87.63044, 41.79625), proj4string=CRS("+proj=longlat +datum=WGS84"))
SP <- spTransform(SP, CRS("+proj=utm +north +zone=16T +datum=WGS84"))
더 의미가있는 경우 이벤트 데이터를 원래 형식으로 작업하려고합니다. 궁극적으로 시카고 도시 경계를 이벤트 데이터와 대조하여 확인할 수있는 다각형으로 바꾸는 데 도움이 필요합니다. 이 같은:
gContains(city, SP)
포인트가 도시에 있는지 확인할 수 있도록 셰이프 파일과 이벤트 데이터를 공통 (그리고 정확한) 참조 프레임으로 가져 오려면 어떻게해야합니까? 렌더링 할 수있는 모든 도움을 주시면 감사하겠습니다.
일반적인 조언은 다음과 rgdal::readOGR같이 셰이프 파일을 읽는 것입니다.
library(rgdal)
city = readOGR(".", "City_Boundary")
이는 city적절한 CRS를 제공 할뿐만 아니라 :
proj4string(city)
[1] "+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"
재 투영하지 않고 변경하려는 경우에도 경고합니다.
proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84")
Warning message:
In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
A new CRS was assigned to an object with an existing CRS:
+proj=tmerc +lat_0=36.66666666666666 +lon_0=-88.33333333333333 +k=0.9999749999999999 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform in package rgdal
이 일 후에 그 의미 city입니다 잘못된 . city를 사용하여 재 투영 합니다 spTransform.
city = spTransform(city, CRS("+proj=utm +north +zone=16T + datum=WGS84"))
귀하의 rgeos::gContainsUTM에 대한 두 CRS 문자열이 다르기 때문에 호출이 실패 할 수 있습니다; + datum=WGS84전자 의 공백을 제거하면 괜찮습니다 ( TRUE).
이 기사는 인터넷에서 수집됩니다. 재 인쇄 할 때 출처를 알려주십시오.
침해가 발생한 경우 연락 주시기 바랍니다[email protected] 삭제
몇 마디 만하겠습니다