Thank you for your rapid response, Edzer. I thought my maths on the coordinates was a bit dodgy and my understanding of the deep depths of crs's is rather superficial, so the explanation is appreciated. I shall dig through your workshop and see where I end up.
kind regards dave ----- Original Message ----- From: "Edzer Pebesma" <edzer.pebe...@uni-muenster.de> To: r-sig-geo@r-project.org Sent: Saturday, February 11, 2023 12:43:15 AM Subject: Re: [R-sig-Geo] Lo Cape projections in sf and terra I ran into this a few months ago while preparing for a workshop I gave (online) to SASA22 participants in South Africa: https://edzer.github.io/SASA22/workshop.html If you divide the geometry by -1 (or multiply it by -1) I don't think you can simply set back the CRS to the old value (EPSG:22289) as that really assumes westing/southing; the relevant part of the WKT2 representation is: CS[Cartesian,2], AXIS["westing (Y)",west, ORDER[1], LENGTHUNIT["metre",1]], AXIS["southing (X)",south, ORDER[2], LENGTHUNIT["metre",1]], It is intentional that the CRS disappears after you do math operations on the coordinates. But indeed, most of the plotting mechanism in R for spatial data simply assume the first axis to be easting, the second to be northing. On 10/02/2023 20:55, dave furniss wrote: > 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 -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo