r - 使用 raster 和 sp 库在 R 中将坐标的 data.frame 从 Lambert93 重新投影到 WGS83

Question

我正在尝试将一些坐标从 Lambert 93 重新投影到 WGS84。我花了足够的时间寻找文档并更好地理解，但我还没有看到解决方案。

我正在寻找可以解释我哪里出错的人。

require(sp)
require(raster)
# >>>>>>>>>>>>>>>>> 
# >>>>>>>>>>>>>>>>> I define the coordinates that I will use later
# >>>>>>>>>>>>>>>>> 

# define coordinates ----
# Lambert 93
# epsg 2154
crs_l93 <-  " +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3
+x_0=700000+y_0=6600000 +ellps=GRS80 +units=m +no_defs"

# WGS84
# epsg 4326
crs_wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# >>>>>>>>>>>>>>>>> 
# >>>>>>>>>>>>>>>>> I define the coordinates i want to project (coord_l93)
# >>>>>>>>>>>>>>>>> and the coordinates goal (what i want them to be projected to : the result)
# >>>>>>>>>>>>>>>>> The coordinates where projected on https://epsg.io/ manually.
# >>>>>>>>>>>>>>>>> 

#  get data ----
  # brut
coord_l93 <-   data.frame("coord_x" = c(839500 , 830500 , 826500 , 826500 ) ,
                        "coord_y" = c(6458500, 6461500, 6467500, 6470500))

# cherché sur https://epsg.io/
coord_wgs84_hoped <- data.frame("coord_x" = c(4.7771833 , 4.6633664 , 4.6139595 , 4.6147419 ) ,
                            "coord_y" = c(45.2116938, 45.2404655, 45.2952272, 45.3222348))

# >>>>>>>>>>>>>>>>> 
# >>>>>>>>>>>>>>>>> I set my data as SpatialPointsDataFrame and i set their projection
# >>>>>>>>>>>>>>>>> 

# define new variable
sp_brut_l93 <- coord_l93
sp_brut_wgs84_hoped <- coord_wgs84_hoped 

# define projection 

sp::coordinates(sp_brut_l93) <- sp_brut_l93
proj4string(sp_brut_l93) <- CRS(crs_l93)
sp::coordinates(sp_brut_wgs84_hoped) <- sp_brut_wgs84_hoped
proj4string(sp_brut_wgs84_hoped) <- CRS(crs_wgs84)


# >>>>>>>>>>>>>>>>> 
# >>>>>>>>>>>>>>>>> I make the projection by two different ways
# >>>>>>>>>>>>>>>>> 

# project sp_coord_l93 into wgs84
sp_proj_wgs84 <- spTransform(spbv, CRS(crs_wgs84))
sp_proj_wgs84_V2 <- spTransform(spbv, "+init=epsg:4326")

# >>>>>>>>>>>>>>>>> 
# >>>>>>>>>>>>>>>>> final check 
# >>>>>>>>>>>>>>>>> 
cat(sp_proj_wgs84)
cat(sp_brut_wgs84_hoped)

# check proj
   # check that the to ways of projecting are identical
identical(coordinates(sp_proj_wgs84), coordinates(sp_proj_wgs84_V2))
   # check that the projected is same as the hoped
identical(coordinates(sp_brut_wgs84_hoped), coordinates(sp_proj_wgs84))

我的预测发生错误，我真的错过了这里的重点。希望任何人都可以在这里解释我。

谢谢你

标签： rgiscartography

Answer 1

这是您的代码，简化且有效

首先用点定义坐标参考系统和data.frames

crs193 <- "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"
wgs84 <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

xyl93 <- data.frame("coord_x" = c(839500 , 830500 , 826500 , 826500 ) ,
                    "coord_y" = c(6458500, 6461500, 6467500, 6470500))

xy84 <- data.frame("coord_x" = c(4.7771833 , 4.6633664 , 4.6139595 , 4.6147419 ) ,
                   "coord_y" = c(45.2116938, 45.2404655, 45.2952272, 45.3222348))

通过识别 data.frame 中保存坐标的列来创建 SpatialPoints（不是通过分配坐标参考系，尽管我们也这样做）

library(sp)
# method 1
coordinates(xyl93) <- ~ coord_x + coord_y
proj4string(xyl93) <- CRS(crs193)

# method 2
xy84 <- sp::SpatialPoints(xy84, proj4string=CRS(wgs84))

现在我们可以改造

x <- spTransform(xy84, crs193)
coordinates(x)
#      coord_x coord_y
#[1,]  839500 6458500
#[2,]  830500 6461500
#[3,]  826500 6467500
#[4,]  826500 6470500