当两个空间对象具有不同的 CRS 时绘制空间数据 [英] Plotting spatial data when two spatial objects have different CRS

问题描述

我有一个空间多边形对象和一个空间点对象.后者是从数据帧中的 xy 纬度数据向量(分别称为纬度和经度)创建的，而前者只是使用 rgdal 直接读入 R.我的代码如下:

I have a spatial polygon object and a spatial points object. The latter was created from xy latlong data vectors (called latitude and longitude, respectively) in a dataframe, while the former was simply read into R directly using rgdal. My code is as follows:

boros <- readOGR(dsn = ".", "nybb")
rats <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv", header = TRUE)
coordinates(rats) <- ~longitude + latitude

此时没有空间对象被投影.如果我按如下方式投影这些对象:

At this point neither spatial object is projected. If I project these objects as follows:

proj4string(boros) <- CRS("+proj=lcc")
proj4string(rats) <- CRS("+proj=lcc")

现在两个对象都被投影了，并且都将成功映射到 plot() 函数，如下所示:

Both objects are now projected, and both will successfully map with the plot() function as follows:

plot(boros)
plot(rats)

但是，当我尝试将它们绘制在一起时:

However when I try to plot them together:

plot(boros)
plot(rats, add = TRUE)

我只得到第一个情节，没有将老鼠对象叠加在 boros 上.然而，这是一个大问题，我没有收到错误消息，所以我一直无法确定这两个能够相互交谈的空间对象之间的断开连接在哪里.这两个命令都可以顺利运行，没有错误或警告，但我只剩下一个情节.当我使用 proj4string() 检查每个对象的投影时，我会为每个对象返回相同的投影.

I get the first plot only, without the rats object superimposed over boros. However, and this is the big problem, I get NO error message, so I have been unable to determine where the disconnection is between these two spatial objects being able to speak to each other. Both commands run smoothly without error or warning, yet I am left with just the single plot. And when I check the projections of each object with proj4string() I get the same projection returned for each object.

我现在已经在几天内花了很多很多小时尝试各种方法来创建两个空间对象，它们的 CRS 和投影匹配，以便可以使用 plot() 将它们映射在一起.顺便说一下，我采用的一种方法是在 ArcGIS 中为老鼠对象创建一个 shapefile，它可以很好地创建文件.但是我仍然无法让两个空间对象在 R 中协同工作.我为这两个对象尝试了许多不同的投影，两个对象上的 spTransform，似乎没有任何效果.非常感激任何的帮助.我还包含了一个 Dropbox 链接，其中包含我上面描述的 2 个数据文件:https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384

I have now spent many, many hours over several days trying various ways of creating two spatial objects whose CRS and projections match such that they can be mapped together using plot(). Incidentally, one approach I took was to create a shapefile in ArcGIS for the rats object, which worked fine to create the file. But I am still left with the same inability of the two spatial objects to work together in R. I have tried many different projections for both objects, spTransform on both objects, nothing seems to work. Any help would be most appreciated. I have also included a dropbox link with the 2 data files I have described above: https://www.dropbox.com/sh/x0urdo6guprnw8y/tQdfzSZ384

推荐答案

因此，正如一些评论指出的那样，问题在于您的数据和地图具有不同的投影.

So, as some of the comments point out, the problem is that your data and your maps have different projections.

看起来您的地图来自 纽约市教育部城市规划.WGS84 (longlat) 中的 shapefile 绝对不是，但是文件中不包含 CRS(顺便说一下，这非常令人失望...).然而，有一个 元数据文件 表明 shapefile预计 EPSG 2263.

It looks like your map comes from the NYC Department of City Planning. The shapefile is definitely not in WGS84 (longlat), but the CRS is not included in the file (which is very disappointing by the way...). Nevertheless, there is a metadata file which indicates that the shapefile is projected as EPSG 2263.

为了在 R 中使用它，我们需要一个投影字符串.在 R 中获得它的惯用方法是:

In order to make use of this in R we need a projection string. The idiomatic way to get this in R is:

library(rgdal)
EPSG <- make_EPSG()
NY   <- with(EPSG,EPSG[grepl("New York",note) & code==2263,]$prj4)
NY
# [1] "+proj=lcc +lat_1=41.03333333333333 +lat_2=40.66666666666666 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000.0000000001 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"

现在我们可以将您的地图重新投影到 WGS84 中，或者将您的数据重新投影到地图 CRS 中.

Now we can either take your map and reproject that into WGS84, or take your data and reproject that into the map CRS.

setwd("< directory with all your files >")
data   <- read.csv("nycrats_missing_latlong_removed_4.2.14.csv")

# First approach: reproject map into long-lat
wgs.84       <- "+proj=longlat +datum=WGS84"
map          <- readOGR(dsn=".",layer="nybb",p4s=NY)
map.wgs84    <- spTransform(map,CRS(wgs.84))
map.wgs84.df <- fortify(map.wgs84)
library(ggplot2)
ggplot(map.wgs84.df, aes(x=long,y=lat,group=group))+
  geom_path()+
  geom_point(data=data, aes(x=longitude,y=latitude, group=NULL),
             colour="red", alpha=0.2, size=1)+
  ggtitle("Projection: WGS84 longlat")+
  coord_fixed()

# Second approach: reproject data
map.df <- fortify(map)
rats <- SpatialPoints(data[,c("longitude","latitude")],proj4string=CRS(wgs.84))
rats <- spTransform(rats,CRS(NY))
rats.df <- data.frame(coordinates(rats))
ggplot(map.df, aes(x=long,y=lat,group=group))+
  geom_path()+
  geom_point(data=rats.df, aes(x=longitude,y=latitude, group=NULL),
             colour="red", alpha=0.2, size=1)+
  ggtitle("Projection: NAD83.2263")+
  coord_fixed()

中央公园没有老鼠?

这篇关于当两个空间对象具有不同的 CRS 时绘制空间数据的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！