::list_coordOps("EPSG:21781", "EPSG:2056") rgdal
When I started working with swiss geodata a couple of years ago, Switzerland was beginning to transition from it’s old coordinate system CH1903 LV03
to CH1903+ LV95
(EPSG 2056).
The new coordinate system was awkward from the start. First off, I suddenly was un unable to convert data from WGS84 using a custom R script that swisstopo had previously provided for CH1903 LV03
. The answer to my question on a swisstopo google group (in 2017) was to use swisstopo’s REST api.
This problem was seemingly solved with the emergence of sf
and it’s built-in methods to reproject coordinates easily. Until one day I realized that the coordinate transformations to and from EPSG 2056 were fairly imprecise. The transformations did not account for the spatially varying offset for which the new coordinate system was implemented in the first place.
I filed an issue on sf’s github repo beginning of 2020, but it wasn’t until recently (with PROJ 7.0.0) that I was able to precisely transform my data to and from EPSG 2056.
The issue is highlighted with the function list_coordOps()
from rgdal
. Apparently, the transformation was lacking a grid with the name ch_swisstopo_CHENyx06a.tif
(the situation is now different, see my edit below).
Candidate coordinate operations found: 2
Strict containment: FALSE
Visualization order: TRUE
Source: EPSG:21781
Target: EPSG:2056
Best instantiable operation has only ballpark accuracy
Description: Inverse of Swiss Oblique Mercator 1903M + Ballpark geographic
offset from CH1903 to CH1903+ + Swiss Oblique Mercator
1995
Definition: +proj=pipeline +step +inv +proj=somerc +lat_0=46.9524055555556
+lon_0=7.43958333333333 +k_0=1 +x_0=600000+y_0=200000 +ellps=bessel
+step +proj=somerc +lat_0=46.9524055555556 +lon_0=7.43958333333333
+k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel
Operation 1 is lacking 1 grid with accuracy 0.2 m
Missing grid: ch_swisstopo_CHENyx06a.tif
URL: https://cdn.proj.org/ch_swisstopo_CHENyx06a.tif
With my current version of sf
and rgdal
, I need to download this grid manually and move it to the PROJ directory. Running list_coordOps
now shows a different output.
::list_coordOps("EPSG:21781", "EPSG:2056") rgdal
Candidate coordinate operations found: 2
Strict containment: FALSE
Visualization order: TRUE
Source: EPSG:21781
Target: EPSG:2056
Best instantiable operation has accuracy: 0.2 m
Description: Inverse of Swiss Oblique Mercator 1903M + CH1903 to
CH1903+ (1) + Swiss Oblique Mercator 1995
Definition: +proj=pipeline +step +inv +proj=somerc
+lat_0=46.9524055555556
+lon_0=7.43958333333333 +k_0=1 +x_0=600000
+y_0=200000 +ellps=bessel +step
+proj=hgridshift
+grids=ch_swisstopo_CHENyx06a.tif +step
+proj=somerc +lat_0=46.9524055555556
+lon_0=7.43958333333333 +k_0=1
+x_0=2600000 +y_0=1200000 +ellps=bessel
The downloaded file is a so called transformation grid: A raster dataset in WGS84 containing information on the lat and lon offset for a given cell.
library(tmap)
library(sf)
library(terra)
library(cowplot)
<- terra::rast("ch_swisstopo_CHENyx06a.tif")
ch_swisstopo_CHENyx06a # Swissboundaries3D data from swisstopo
<- read_sf("data-git-lfs/swissboundaries/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp")
landesgebiet
<- st_transform(landesgebiet, 2056)
landesgebiet
tm_shape(ch_swisstopo_CHENyx06a) +
tm_raster(style = "cont",palette = "RdYlBu") +
tm_shape(landesgebiet) + tm_borders()
I now can use this grid to precisely transform coordinates to and from ESPG 2056. I will demonstrate this by visualizing the directional offset, similar to the transformation grid I downloaded.
library(dplyr)
library(ggplot2)
library(tidyr)
<- function(h){
heading_north case_when(
>=0 & h<=90~90-h,
h>90~450-h,
hTRUE~abs(h) + 90
)
}
<- function(cellsize, land){
make_grid
<- st_transform(land, 2056)
land
<- st_make_grid(land, cellsize, what = "centers")
gri
<- st_transform(gri, 21781)
gri2
<- cbind(
mat st_coordinates(gri),
st_coordinates(gri2)
)
colnames(mat) <- c("E", "N", "X", "Y")
%>%
mat ::as_tibble() %>%
tibblemutate(
E_norm = E - (X+2e6),
N_norm = N - (Y+1e6),
dist = (E_norm^2+N_norm^2)^0.5,
heading = atan2(N_norm, E_norm),
h = heading*180/pi,
heading_deg = heading_north(h)
)
}
# Add some cities as labels
<- tribble(
markers ~E, ~N, ~name, ~nudge_x, ~nudge_y,
2598633.750, 1200386.750, "Bern",15000,15000,
2682217.000, 1247945.250 , "Zurich",15000,-15000,
2598410, 1185268, "apparent epicenter",15000,-15000,
)
# Make my custom grid, similar to the transformation grid
<- make_grid(1000, landesgebiet)
df
# https://coolors.co/ef476f-ffd166-06d6a0-118ab2-073b4c
<- paste0("#",c("EF476F","FFD166","06D6A0", "118AB2","EF476F")) mycols
<- make_grid(4000, landesgebiet) %>%
mygrid st_as_sf(coords = c("E","N"), remove = FALSE) %>%
st_set_crs(2056)
<- mygrid[landesgebiet, , op = st_intersects]
mygrid
<- Sys.getenv("plot_bg_col")
plot_bg_col <- Sys.getenv("text_col")
text_col
<- ggplot(mygrid, aes(E, N, color = heading_deg)) +
main_plot # geom_sf(data = bb_ch, inherit.aes = FALSE, fill = plot_bg_col, color = "NA") +
geom_point(size = .4) +
geom_spoke(aes(angle = heading, radius = scales::rescale(dist, 5e3, 10e3))) +
scale_color_gradientn(colors = mycols) +
labs(title = "Warping Switzerland back into shape",
subtitle = "How off were the old coordinates?",
caption = "The spokes show the offset between EPSG 21781 and EPSG 2056.\nData from swisstopo. Visualized by Nils Ratnaweera") +
coord_equal() +
theme_void() +
theme(legend.position = "none",
plot.title = element_text(hjust = .5),
plot.subtitle = element_text(hjust = .5),
text = element_text(color = text_col),
plot.margin = margin(10,10,10,10),
panel.background = element_rect(colour = NA, fill = plot_bg_col),
plot.background = element_rect(fill = plot_bg_col,color = NA),aspect.ratio = 2/3)
<- function(deg){deg*pi/180}
deg2rad
<- tibble(
legend_df angle_deg1 = seq(0,359,10),
angle_deg2 = heading_north(angle_deg1),
angle_rad = deg2rad(angle_deg2),
x = cos(angle_rad),
y = sin(angle_rad)
)
= 2
fac <- tibble(
windrad angle_deg = c(0,90,180,270),
text = c("N","E","S","W"),
angle_rad = deg2rad(heading_north(angle_deg)),
x = cos(angle_rad)*fac,
y = sin(angle_rad)*fac
)
<- ggplot(legend_df, aes(x, y, color = angle_deg1)) +
legend_plot geom_point(size = .8) +
geom_text(data = windrad, aes(label = text), color = text_col, size = 2) +
geom_spoke(aes(angle = angle_rad,radius = .6)) +
scale_color_gradientn(colors = mycols) +
coord_equal() +
theme_void() +
theme(legend.position = "none")
ggdraw(main_plot) +
draw_plot(legend_plot, .95, .95, .2, .2,
hjust = 1,vjust = 1, halign = 1, valign = 3) +
theme(panel.background = element_rect(fill = plot_bg_col))
ggsave("preview.jpg", width = 15, height = 10, units = "cm")
EDIT (2022-04-24):
At the time of writing this blogpost, the bash command to download and move the file to the correct place was as follows:
wget https://cdn.proj.org/ch_swisstopo_CHENyx06a.tif`
sudo mv ch_swisstopo_CHENyx06a.tif /usr/share/proj/
Today however, I needed to download a zip file and unzip it to /usr/share/proj/
. I assume this will change again in the near future.
https://download.osgeo.org/proj/proj-datumgrid-europe-1.5.zip
sudo unzip proj-datumgrid-europe-1.5 /usr/share/proj/