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

Я новичок в ГИС в R, и я пытаюсь добавить координаты широты и долготы или координаты 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")

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

Моя конечная цель - определить, происходило ли в Чикаго несколько событий с геотегами. Мне удобно конвертировать те точки, которые указаны по широте / долготе, в 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)

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


r gis rgdal maptools
person user2047457    schedule 25.02.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Общий совет - читать шейп-файлы с 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::gContains может завершиться неудачно, потому что две строки CRS для utm различаются; удалите пробел в + datum=WGS84 в первом, и все будет в порядке (TRUE).

person Edzer Pebesma    schedule 25.02.2015
comment
Спасибо, отличный ответ! Именно то, что мне нужно. Я использовал proj4string from city, чтобы преобразовать мои точки широты и долготы в тот же UTM, что и город. СП ‹- spTransform (СП, city @ proj4string) - person user2047457; 25.02.2015