Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Created June 15, 2025 21:20
Show Gist options
  • Save mdsumner/e5cf8d10574bf61f99c6a84f5dc1378e to your computer and use it in GitHub Desktop.
Save mdsumner/e5cf8d10574bf61f99c6a84f5dc1378e to your computer and use it in GitHub Desktop.
## we can get that WKT by using Oblique Mercator params 
prj <- "+proj=omerc +lonc=-122.6 +lat_0=40 +gamma=21 +x_0=1000000"

shp <- "/vsicurl/https://github.com/noaa-nwfsc/VMS-pipeline/raw/refs/heads/main/spatial_data/map_rotation/fivekm_grid_extent_rect_21.shp"
library(terra)
#> terra 1.8.52
v <- vect(shp)

wkt0 <- crs(prj)

## round trip with out proposed PROJ string form of the wkt
plot(project(project(v, "EPSG:4326"), wkt0))
plot(v, add = T, lwd = 2)  ## it's the same

## get some raster data
dsn <- "/vsicurl/https://data.source.coop/alexgleith/gebco-2024/GEBCO_2024.tif"
target <- rast(ext(v), res = 2500, crs = prj)  ## tweak resolution to be sensible (not too many pixels)
target
#> class       : SpatRaster 
#> size        : 990, 826, 1  (nrow, ncol, nlyr)
#> resolution  : 2499.146, 2500.01  (x, y)
#> extent      : -231257, 1833038, -1252242, 1222769  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=omerc +lat_0=40 +lonc=-122.6 +alpha=0 +gamma=21 +k=1 +x_0=1000000 +y_0=0 +datum=WGS84 +units=m +no_defs

topo <- project(rast(dsn), target, by_util = TRUE)
plot(topo); plot(v, add = TRUE)

Created on 2025-06-15 with reprex v2.0.2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment