library(sf)
## Linking to GEOS 3.5.0, GDAL 2.1.1, proj.4 4.9.3
prj <- "+proj=omerc +datum=WGS84 +lonc=165 +lat_0=-22 +alpha=23 +gamma=23"
crs <- st_crs(prj)
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
data(wrld_simpl)
w <- st_transform(st_as_sf(subset(wrld_simpl, NAME == "Australia")), crs)
plot(st_geometry(w))
plot(st_geometry(st_graticule(st_geometry(w))), add = T)
## sampling problems
plot(st_geometry(st_graticule(st_bbox(w) + c(-1, -1, 1, 1) * 5e6, crs = st_crs(w))))
## cannot find an option to get this to work (is there a better approach?)
#plot(st_geometry(st_graticule(st_bbox(w) + c(-1, -1, 1, 1) * 1e7, crs = st_crs(w))))
#Error in apply(dxdy, 1, function(x) atan2(x[2], x[1]) * 180/pi) :
# dim(X) must have a positive length
## try very early "finite-element" method, using contour()
#devtools::install_github("mdsumner/graticule@auto-graticule")
library(graticule)
plot(graticule:::auto_graticule(as(w, "Spatial")))
ww <- as(raster::extent(as(w, "Spatial")) + 5e6, "SpatialPolygons")
proj4string(ww) <- CRS(prj)
plot(graticule:::auto_graticule(as(ww, "Spatial")))
w1 <- spTransform(wrld_simpl, CRS(prj))
## works ok for the whole world
plot(graticule:::auto_graticule(w1))
## but map-wise we need to decompose to primitives, trim out wrapping elements
## that are more complicated than the usual -180/180 problem
plot(graticule:::auto_graticule(w1), col = "darkgrey")
plot(w1, add = TRUE, border = "firebrick")
## see the over-sampling at 180
## need to switch coords for contouring and then fakey back to real
w2 <- spTransform(subset(wrld_simpl, NAME %in% c("Australia", "Antarctica")), CRS(prj))
plot(graticule:::auto_graticule(w2))
plot(w2, add = TRUE, border = "dodgerblue")
