CODEWITHSUNDEEP
rgiscartography
Lilian Villenave
Lilian Villenave

Reputation: 1

Reproject data.frame of coordinates from Lambert93 to WGS83, in R, using raster, and sp library

I am trying to reproject some coordinates from Lambert 93 to WGS84. I spent enought time looking for documentation and to understand better but I don't see a solution yet.

I'm looking someone that could explain me where I got wrong.

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

My projection happen wrongly and I really miss the point here. Hope anyone can explain me here.

Upvotes: 0

Views: 1170

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47146

Here is your code, simplified, and working

First define coordinate reference systems and data.frames with points

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

Create SpatialPoints by identifying the columns in the data.frame that hold the coordinates (not by assigning the coordinate reference system, although we also do that)

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

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

Now we can transform

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

Upvotes: 1

Related Questions