Hello all
I'm struggling with a projection system used by miners in South Africa, the Lo
Cape system. I've put together some commented dummy code that I hope
illustrates the issues. If someone could give me some guidance on how to
convert between this and other crs's, that would be much appreciated.
### I have lidar collected in an old South African coordinate system,
### the Lo Cape system, which I believe is EPSG 22289 for my site. I've made
### some dummy data to illustrate the problems I'm having in converting
### data to this projection system or from this to another crs.
library(sf)
library(terra)
### Make a point from GoogleEarth coordinates and set projection to wgs84
### and then convert to Lo 29 Cape projection (EPSG 22289)
(bridge <- st_sf(st_sfc(st_point(c(29.472945, -25.982804)), crs = 4326)))
(bridgest <- st_transform(bridge, crs = 22289))
### Notice the coordinates have their signs reversed. So the point is
### now northern and western hemisphere when it should be in the
### southern and eastern hemisphere. For a set of points I can
### correct this by dividing the goemetry by -1 as below
(bridgeT <- st_set_geometry(bridgest, st_geometry(bridgest)/-1))
### However, the crs is now lost but the signs on the
### coordinates are switched to what they should be
st_crs(bridgeT) <- 22289
bridgeT
### Lidar data was collected in the Lo 29 cape dataum. I've made
### a dummy raster of this data from the original data:
dumrast <- rast(ncol = 7, nrow = 6, xmin = 47382, xmax = 47389,
ymin = -2874721, ymax = -2874715, crs = "EPSG:22289")
values(dumrast) <- c(1562.78, 1562.58, 1562.50, 1562.99, 1563.27, 1563.81,
1564.33, 1562.81,
1562.66, 1562.61, 1563.35, 1563.57, 1563.96, 1564.38, 1562.87, 1562.72,
1563.57, 1563.59, 1563.78, 1564.01, 1564.27, 1562.92, 1562.81, 1563.85,
1563.76, 1563.83, 1564.26, 1564.24, 1562.73, 1562.62, 1562.89, 1563.11,
1563.21, 1563.47, 1564.38, 1562.69, 1562.55, 1562.47, 1562.80, 1563.05,
1563.76, 1563.98)
dumrast
### If we plot the DEM and the point, they match
plot(dumrast)
plot(bridgeT, add = TRUE, col = 'red')
### If I try and convert the point crs from Lo Cape to wgs 84,
### the negative is lost and now instead of being in the southern
### and eastern hemisphere, the point is in the northern and eastern
### hemisphere
bridgeT
(pnt <- st_transform(bridgeT, 4326))
### I derived a drainage from the DEM in ArcGIS. ArcGIS doesn't seem to
### recognize the crs when the tif is saved in r. Nonetheless I've used
### a portion of the drainage to create a line segment which should match
### another known point built from Google Earth coordinates:
(seep <- st_sf(st_sfc(st_point(c(29.459556, -26.004615)), crs = 4326)))
### As before, converting the crs from wgs 84 to Lo Cape again switches
### the positive and negative values
(seept <- st_transform(seep, 22289))
(seepT <- st_set_geometry(seept, st_geometry(seept)/-1))
st_crs(seepT) <- 22289
seepT
### Now I make a short section of the drainage line derived from the
### lidar DEM
m <- matrix(c(46056.5,-2877155, 46048.5,-2877146, 46046.5,-2877142,
46043.5,-2877140, 46037.5,-2877133, 46037.5,-2877131, 46042.5,-2877126,
46045.5,-2877118, 46054.5,-2877108, 46057.5,-2877099, 46060.5,-2877097),
ncol=2, byrow = TRUE)
(ditch <- st_sf(st_sfc(st_linestring(m)), crs = 22289))
plot(ditch, axes = TRUE)
plot(seepT, add = TRUE, col = 'red')
### They match close enough for now.
### If I try and transform the linestring to wgs coordinates,
### I again lose the negative in the coordinates
(ditch_wgs <- st_transform(ditch, 4326))
### reprojecting a dem
(demWgs <- project(dumrast, "epsg:4326"))
### The crs for the raster also loses its negative values and is rotated.
### So please can someone provide guidance on how best to do these
### transformations or suggest references that may help me find better
### solutions than my hacks.
> sessionInfo()
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=English_South Africa.utf8 LC_CTYPE=English_South Africa.utf8
[3] LC_MONETARY=English_South Africa.utf8 LC_NUMERIC=C
[5] LC_TIME=English_South Africa.utf8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] terra_1.6-47 sf_1.0-9
loaded via a namespace (and not attached):
[1] Rcpp_1.0.9 magrittr_2.0.3 units_0.8-1 tidyselect_1.2.0
[5] R6_2.5.1 rlang_1.0.6 fansi_1.0.3 dplyr_1.0.10
[9] tools_4.2.2 grid_4.2.2 KernSmooth_2.23-20 utf8_1.2.2
[13] cli_3.5.0 e1071_1.7-12 DBI_1.1.3 class_7.3-20
[17] assertthat_0.2.1 tibble_3.1.8 lifecycle_1.0.3 codetools_0.2-18
[21] vctrs_0.5.1 glue_1.6.2 proxy_0.4-27 compiler_4.2.2
[25] pillar_1.8.1 generics_0.1.3 classInt_0.4-8 pkgconfig_2.0.3
>
Thank you
Dr David Furniss, Wits University.
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
