Impossibility to transform coordinates in degree with sf [R]

Question

I would like to convert coordinates from two different CRS on R with the package sf

	From	To
DATUM	ETRS89	BELGIAN DATUM 72
Ellipsoid	GRS80	HAYFORD24
Projection	None	LAMBERT [1972]
Coordinates	GEOGRAPHIC	PLANE
Unit	degree+minute	meter

The conversion from data with degree to meter works properly with this code :

d = st_as_sf(d, coords=c("lon","lat"), crs="EPSG:4258" )
dt = st_transform(d, "EPSG:31370")

However, when I use it with degree + minutes, I obtain empty points...

(As example, a point with lat: 536.1314 and long: 4943.265 should give an output of x: 229332.264 y: 525235.654)

The problematic dataset:

d <- structure(list(lat = c(4943.26460468966, 4943.264349, 4943.264298, 
4943.26419, 4943.26397666667, 4943.2637825), lon = c(536.131357103448, 
536.1306775, 536.13063, 536.13047, 536.130253333333, 536.1299975
), z = c(400.473048275862, 400.2688, 400.183, 400.08975, 399.863333333333, 
399.82925)), row.names = c(NA, 6L), class = "data.frame")

I am also joining an example of data also in ETRS89 but which is working properly (in degree):

d <- structure(list(name = c("12112", "12201", "12202", "12203", "12208", 
"12209"), lat = c(50.0855877916667, 50.0855997333333, 50.0856167121212, 
50.0856117666667, 50.0856332380952, 50.0855853833333), lon = c(4.55540004166667, 
4.55512903333333, 4.55519089393939, 4.5551871, 4.55524742857143, 
4.55543285), id = c("12112", "12201", "12202", "12203", "12208", 
"12209")), row.names = c(NA, 6L), class = "data.frame")

I do not know how to signify the sf::st_as_sf(), I am in degree + minutes

EDIT

As example, a point with lat: 536.1314 and long: 4943.265 should give an output of x: 238942.84 y: 46236.164 (using the software cConvert)

Answer 1

The numbers in your data frame are not in degrees:

> d
       lat      lon        z
1 4943.265 536.1314 400.4730
2 4943.264 536.1307 400.2688

It looks like 4943.265 should be 49 degrees 43.265 minutes. You'll need to convert this to decimal degrees which will be 49 + (43.265/60) which is 49.72108.

To compute this you'll need to take the residue and modulus by 100 and work it out. Here's a useful function:

getdec = function(x){x%/%100 + (x%%100)/60}

then do:

d$latfix = getdec(d$lat)
d$lonfix = getdec(d$lon)

Now make a spatial data frame

s = st_as_sf(d, coords=c("lonfix","latfix"), crs=4258)

and transform:

> st_transform(s, 31370)
Simple feature collection with 6 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 238941.1 ymin: 46233.94 xmax: 238942.7 ymax: 46235.49
Projected CRS: Belge 1972 / Belgian Lambert 72
       lat      lon        z                  geometry
1 4943.265 536.1314 400.4730 POINT (238942.7 46235.49)
2 4943.264 536.1307 400.2688    POINT (238941.9 46235)
3 4943.264 536.1306 400.1830 POINT (238941.9 46234.91)
4 4943.264 536.1305 400.0897 POINT (238941.7 46234.71)
5 4943.264 536.1303 399.8633 POINT (238941.4 46234.31)
6 4943.264 536.1300 399.8293 POINT (238941.1 46233.94)

Which is close enough to what you expected.

I'm not sure this will work for negative numbers, but you should be okay with Belgium which is conveniently all in the (+,+) quadrant.