3

This is the data I have in input:

> points
# A tibble: 4 × 2
  latitude   longitude
  <chr>      <chr>    
1 49.5150946 3.4798894
2 49.6635237 3.2990396
3 49.1200986 3.3175965
4 49.5638522 3.6248577

After some research I have come with this attempt that does not work:

city_proj <- "+proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
points <- st_as_sf(points, coords = c(1, 2))
st_crs(points) <- 2154
(new_points <- st_transform(points, crs = city_proj))

it leads to:

Simple feature collection with 4 features and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -1939867 ymin: -4828373 xmax: -1939867 ymax: -4828373
CRS:           +proj=lcc +lat_1=43 +lat_2=62 +lat_0=30 +lon_0=10 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs
# A tibble: 4 × 1
             geometry
*         <POINT [m]>
1 (-1939867 -4828373)
2 (-1939867 -4828373)
3 (-1939867 -4828373)
4 (-1939867 -4828373)

I can see that the CRS is not okay, but I do not know how to set it up.

2
  • 1
    The proj4 definition I can find for Lambert-93 is different from the one you are using: '+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs' Source
    – Kasper
    Dec 22, 2022 at 16:09
  • The CRS you are using is one for all of Europe--3034, I believe.
    – mkennedy
    Jan 17, 2023 at 4:23

1 Answer 1

2

I finally understand how things work together.

Need to know which CRS (coordinate reference system) is used by both sides

So here is the code needed to convert raw latitude/longitude numbers with the sf package:

library(sf)
library(tidyverse)

## raw_coordinates as numbers from geocoding Google maps API
(
  tribble(
    ~ latitude,  ~ longitude,
    49.5150946,    3.4798894,
    49.6635237,    3.2990396,
    49.1200986,    3.3175965,
    49.5638522,    3.6248577
) -> raw_coordinates)

(  
  raw_coordinates
  %>% mutate(
    ## First transform them in geometric points
    st_point = map2(raw_coordinates$longitude,
                    raw_coordinates$latitude,
                    ~ c(.x, .y))
               %>% map(st_point),
    ## Create feature telling we are using ESPG:4236
    google_map_features = st_sfc(st_point, crs = 4236),
    ## Convert to ESPG:2154
    lambert93_features  = st_transform(google_map_features, crs = 2154)
  )
  %>% pull(lambert93_features)
)


which produces:

Geometry set for 4 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 721601.7 ymin: 6891190 xmax: 745225.1 ymax: 6951643
Projected CRS: RGF93 v1 / Lambert-93
POINT (734765.9 6935191)
POINT (721601.7 6951643)
POINT (723185.3 6891190)
POINT (745225.1 6940689)
2
  • if the the data returned is long/lat then its not EPSG:4326. EPSG:4326 is lat/long
    – nmtoken
    Dec 23, 2022 at 14:35
  • I have changed that in the text
    – pietrodito
    Dec 23, 2022 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.