[R-sig-Geo] Lo Cape projections in sf and terra

dave furniss on||ne445934 @end|ng |rom te|kom@@@net
Sat Feb 11 09:22:08 CET 2023 

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.pebesma using uni-muenster.de>
To: r-sig-geo using 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 using 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 using r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

More information about the R-sig-Geo mailing list