R 将 CSV 映射到 SHP
R Map CSV & SHP
我正在尝试在 shapefile 上绘制纬度和经度坐标。我被告知它们具有相同的坐标系 (NAD27)。在研究它之后,我尝试了很多方法,但没有一种是完美的。他们要么映射一个,要么映射另一个,但不能同时映射两者。下面是读入它们的示例,但我什至不知道如何开始绘图,因为我的 csv 来自 sp 包,我的 shp 来自 maptools 包。
pts <- read.csv("locations.csv")
pts <- pts[,4:5]
names(pts) <- c("Latitude", "Longitude")
pts <- pts[complete.cases(pts),]
pts <- SpatialPoints(pts)
border <- maptools::readShapePoly("landman_one shape.shp")
如有任何帮助,我们将不胜感激。
下面是不同的结构。
> str(pts)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
..@ coords : num [1:372, 1:2] 30.9 30.9 31 31 31 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:372] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "Latitude" "Longitude"
..@ bbox : num [1:2, 1:2] 29.5 -105.2 36 -76.1
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "Latitude" "Longitude"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA
> str(border)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 1 obs. of 2 variables:
.. ..$ Id : int 0
.. ..$ PERIM_GEO: num 998432
.. ..- attr(*, "data_types")= chr [1:2] "N" "F"
..@ polygons :List of 1
.. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. .. .. ..@ Polygons :List of 1
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 922511 545445
.. .. .. .. .. .. ..@ area : num 3.45e+10
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:252, 1:2] 954182 954073 954006 953914 953828 ...
.. .. .. ..@ plotOrder: int 1
.. .. .. ..@ labpt : num [1:2] 922511 545445
.. .. .. ..@ ID : chr "0"
.. .. .. ..@ area : num 3.45e+10
..@ plotOrder : int 1
..@ bbox : num [1:2, 1:2] 819851 386743 1042230 669865
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA
另外,我对投影不太熟悉,所以如果有任何指导,我们将不胜感激。另外,如果有办法用同一个包做到这一点,我会接受的。
问题是您的两个空间数据集都没有投影集并且使用的是不同的坐标系。您可以在 @proj4string
行中看到这一点,它们目前都是 NA
。确定坐标系不同的快速方法是使用 @bbox
插槽 ('bounding box')。在 pts
中,它们是 longitude/latitude,因为它们分别在 -180/+180 和 -90/+90 范围内。另一方面,border
具有完全不同单位的值 (819851 386743 1042230 669865
)。这些看起来都不在 NAD27
坐标系中。
为了解决这个问题,首先我会使用 rgdal::readOGR()
或 sf::read_sf()
来读取您的 border
形状文件。如果存在的话,这些函数会加载到一个投影中,所以这应该会告诉您这些数据在哪个投影中。
其次,要为 pts
设置投影,您可以使用 proj4string()
:
proj4string(pts) = CRS("+init=epsg:4326")
一旦 pts
和 border
都有投影集,您可以使用 spTransform()
将一个坐标系转换到另一个坐标系。因为我不知道 border
的 CRS/projection 是什么,所以我会将其转换为 WGS84(pts
的投影),但只需交换它们即可将 pts
转换为与 border
相同的投影,一旦你知道它的 CRS。
border = spTransform(border, CRS("+init=epsg:4326"))
一旦在同一个投影中,它们应该一起绘制。您可以确认他们的预测匹配:
proj4string(border) == proj4string(pts)
我正在尝试在 shapefile 上绘制纬度和经度坐标。我被告知它们具有相同的坐标系 (NAD27)。在研究它之后,我尝试了很多方法,但没有一种是完美的。他们要么映射一个,要么映射另一个,但不能同时映射两者。下面是读入它们的示例,但我什至不知道如何开始绘图,因为我的 csv 来自 sp 包,我的 shp 来自 maptools 包。
pts <- read.csv("locations.csv")
pts <- pts[,4:5]
names(pts) <- c("Latitude", "Longitude")
pts <- pts[complete.cases(pts),]
pts <- SpatialPoints(pts)
border <- maptools::readShapePoly("landman_one shape.shp")
如有任何帮助,我们将不胜感激。
下面是不同的结构。
> str(pts)
Formal class 'SpatialPoints' [package "sp"] with 3 slots
..@ coords : num [1:372, 1:2] 30.9 30.9 31 31 31 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:372] "1" "2" "3" "4" ...
.. .. ..$ : chr [1:2] "Latitude" "Longitude"
..@ bbox : num [1:2, 1:2] 29.5 -105.2 36 -76.1
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "Latitude" "Longitude"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA
> str(border)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 1 obs. of 2 variables:
.. ..$ Id : int 0
.. ..$ PERIM_GEO: num 998432
.. ..- attr(*, "data_types")= chr [1:2] "N" "F"
..@ polygons :List of 1
.. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. .. .. ..@ Polygons :List of 1
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 922511 545445
.. .. .. .. .. .. ..@ area : num 3.45e+10
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:252, 1:2] 954182 954073 954006 953914 953828 ...
.. .. .. ..@ plotOrder: int 1
.. .. .. ..@ labpt : num [1:2] 922511 545445
.. .. .. ..@ ID : chr "0"
.. .. .. ..@ area : num 3.45e+10
..@ plotOrder : int 1
..@ bbox : num [1:2, 1:2] 819851 386743 1042230 669865
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
.. .. ..@ projargs: chr NA
另外,我对投影不太熟悉,所以如果有任何指导,我们将不胜感激。另外,如果有办法用同一个包做到这一点,我会接受的。
问题是您的两个空间数据集都没有投影集并且使用的是不同的坐标系。您可以在 @proj4string
行中看到这一点,它们目前都是 NA
。确定坐标系不同的快速方法是使用 @bbox
插槽 ('bounding box')。在 pts
中,它们是 longitude/latitude,因为它们分别在 -180/+180 和 -90/+90 范围内。另一方面,border
具有完全不同单位的值 (819851 386743 1042230 669865
)。这些看起来都不在 NAD27
坐标系中。
为了解决这个问题,首先我会使用 rgdal::readOGR()
或 sf::read_sf()
来读取您的 border
形状文件。如果存在的话,这些函数会加载到一个投影中,所以这应该会告诉您这些数据在哪个投影中。
其次,要为 pts
设置投影,您可以使用 proj4string()
:
proj4string(pts) = CRS("+init=epsg:4326")
一旦 pts
和 border
都有投影集,您可以使用 spTransform()
将一个坐标系转换到另一个坐标系。因为我不知道 border
的 CRS/projection 是什么,所以我会将其转换为 WGS84(pts
的投影),但只需交换它们即可将 pts
转换为与 border
相同的投影,一旦你知道它的 CRS。
border = spTransform(border, CRS("+init=epsg:4326"))
一旦在同一个投影中,它们应该一起绘制。您可以确认他们的预测匹配:
proj4string(border) == proj4string(pts)