Добавить координаты lat/long или URM в шейп-файлы в R

Question

Я новичок в GIS в R, и я пытаюсь добавить координаты lat/long или UTM в шейп-файл. Я скачал границы города Чикаго (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 & пытался добавить код UTM для зоны 16T:

city<- readShapeSpatial("C:/Users/Luke.Gerdes/Documents/GIS Files/Chicago/City_Boundary.shp") 
proj4string(city) <- CRS("+proj=utm +north +zone=16T + datum=WGS84") 

Однако результирующие данные для меня не имеют смысла. Когда я смотрю на слоты «коорды» внутри «города», значениями координат являются, например, X = 1092925 и Y = 1944820. Я использовал внешний инструмент (http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html), чтобы найти GPS-координаты. Результаты были: lat = 17,511, long = -81,42. Это где-то между Ямайкой и побережьем Гондураса. Мне кажется, что координаты шейп-файла существуют в их собственной вселенной; как, возможно, имеет название 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) 

Как я могу получить мои шейп файлов и данных событий в общую (и точный) системе отсчета, так что я могу проверить, является ли точки в городе? Любая помощь, которую вы можете оказать, будет с благодарностью.

Answer 1

Общий совет читать шейп с 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 

что означает afte r это делает, city есть не правильный. Вы перепроецировать city с помощью spTransform, как в

city = spTransform(city, CRS("+proj=utm +north +zone=16T + datum=WGS84")) 

Ваш rgeos::gContains вызов может потерпеть неудачу, потому что две CRS строки для UTM различаются; удалите место в + datum=WGS84 в первом, и все будет хорошо (TRUE).