Перепроектируйте data.frame координат из Lambert93 в WGS83, в R, используя растровые и sp библиотеки - PullRequest
Новая
       67

Перепроектируйте data.frame координат из Lambert93 в WGS83, в R, используя растровые и sp библиотеки

0 голосов
/ 01 апреля 2020

Я пытаюсь перепроектировать некоторые координаты от Ламберта 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))

Мои прогнозы случаются неправильно, и я действительно скучаю по сути. Надеюсь, что кто-нибудь может объяснить мне здесь.

Спасибо

1 Ответ

0 голосов
/ 01 апреля 2020

Вот ваш код, упрощенный и работающий

Сначала определите системы координат координат и 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))

Создайте SpatialPoints, идентифицируя столбцы в data.frame, который содержит координаты (не назначая систему координат, хотя мы также делаем это) 

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
...