Географический CRS для несоответствующих данных в R

извините, что беспокою вас по этому поводу, но, может быть, уже 5 часов я схожу с ума по этому вопросу и не могу разобраться.

У меня есть набор данных около 37 000 записей. Каждый из них имеет свои значения долготы и широты. Проверяя общие значения, они варьируются соответственно следующим образом: Широта (-54,4871,70,66344) и Долгота (-177,375, 178,4419). Это абсолютно разумно.

Я создал шейп-файл с этими 37 тысячами точек с помощью ArcGIS: все работает нормально.

Затем мне нужно обработать эти данные с помощью R, команда, которую я использовал для своего кода (пакет maptools):

cells <- readShapeSpatial('RES',IDvar="id_obj", 
                          proj4string=CRS("+proj=longlat +datum=WGS84"))

но R выдает ошибку:

Ошибка в validMethod(as(object, superClass)) : Географический CRS для несоответствующих данных: 2,76663393422e+145

(Я понятия не имею, откуда взялось это число, оно не является частью моего набора данных...)

Читая другие сообщения в этом блоге, кажется, что причиной должны быть недопустимые данные для долготы или широты, но, как я уже упоминал выше, это не относится к моему набору данных.

Пробовал создавать разные шейп-файлы, первый не проецировался, используя несколько проекций (WGS84 Mercator, web mercator...), но ошибка всегда одна и та же...

Спасибо за вашу помощь.


r maptools
person Nemesi    schedule 29.01.2014    source источник
comment
Можете ли вы загрузить свой шейп-файл куда-нибудь (Dropbox ??) и опубликовать ссылку? В противном случае это просто догадки. Я видел эту проблему, когда широта и долгота меняются местами.   -  person jlhoward    schedule 30.01.2014
comment
Большое спасибо за это: вот ссылка dropbox.com/sh/zco8y93ttnnt5dz/ BmRTimCPjV?n=212111396   -  person Nemesi    schedule 30.01.2014

Ответы (1)


arrow_upward
5
arrow_downward

Суть в том, что ваш шейп-файл кажется поврежденным.

Шейп-файл точек имеет два основных раздела: раздел coords с координатами точек и раздел данных с данными «атрибутов» (информация о точках, например регион и страна в вашем случае). В вашем шейп-файле также есть долгота и широта в разделе данных, но они не совпадают:

library(rgdal)
setwd("<directory with shapefile...>")
map <- readOGR(dsn=".", layer="test")

range(map@data$Lat)
# [1] -54.48708  70.66344

range(map@coords[,2])
# [1]  -5.448708e+01  2.766634e+145

Повторное проецирование предполагает преобразование информации в разделе coords, поэтому оно не удалось.

Вот обходной путь, но взламывать SpatialPointsDataFrame не очень хорошая идея:

map@coords <- as.matrix(map@data[c("Lon","Lat")])
map@bbox   <- rbind(range(map@coords[,1]),range(map@coords[,2]))
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
proj4string(map) <- CRS(wgs.84)

library(ggplot2)
gg <- data.frame(map@coords)
ggplot(gg) + 
  geom_point(aes(x=Lon,y=Lat), size=1, alpha=0.5, colour="blue") + 
  coord_fixed()


mercator <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
map.mercator <- spTransform(map,CRS=CRS(mercator))
gg <- data.frame(map.mercator@coords)
ggplot(gg) + 
  geom_point(aes(x=Lon,y=Lat), size=1, alpha=0.5, colour="green") + 
  coord_fixed()

Я бы порекомендовал вам заново создать шейп-файл и повторить попытку.

person jlhoward    schedule 29.01.2014
comment
Большое спасибо. Я действительно ценю. - person Nemesi; 30.01.2014