filledcontour.R

library(oce)

## Warning: package 'oce' was built under R version 3.4.4

## Loading required package: gsw

## Warning: package 'gsw' was built under R version 3.4.4

## Loading required package: testthat

data(coastlineWorld)
data(topoWorld)

## 1. topography
op <- par(mfrow=c(2, 1), mar=c(2, 2, 1, 1))
lonlim <- c(-70, -50)
latlim <- c(40, 50)
topo <- decimate(topoWorld, by=2) # coarse to illustrate filled contours
topo <- subset(topo, latlim[1] < latitude & latitude < latlim[2])
topo <- subset(topo, lonlim[1] < longitude & longitude < lonlim[2])
mapPlot(coastlineWorld, type='l',
        longitudelim=lonlim, latitudelim=latlim,
        projection="+proj=lcc +lat_1=40 +lat_2=50 +lon_0=-60")
breaks <- seq(-5000, 1000, 500)
mapImage(topo, col=oce.colorsGebco, breaks=breaks)


par(op)

## brute force raster
library(raster)

## Loading required package: sp

rtopo <- raster(list(x = topo@data$longitude, y = topo@data$latitude, z = topo@data$z), 
                crs = "+init=epsg:4326")
plot(projectRaster(rtopo, crs = "+proj=lcc +lat_1=40 +lat_2=50 +lon_0=-60"), 
     col=oce.colorsGebco(length(breaks)+ 1), breaks=breaks)

## without loss, convert all pixels to polygons (wasteful since 5 explicit coords every time)
scl <- function(x) (x - min(x, na.rm = TRUE))/diff(range(x, na.rm = TRUE))
## wasteful, creating a colour for every pixels
cols <- oce.colorsGebco(length(breaks)+ 1)[scl(values(rtopo)) * length(breaks) + 1]
plot(sf::st_transform(spex::polygonize(rtopo), "+proj=lcc +lat_1=40 +lat_2=50 +lon_0=-60"), 
     col = cols, border = NA)

## draw from the mesh incrementally, so we never make more than the mesh set
## and each drawn polygon uses only 4 from that shared pool
qtopo <- quadmesh::quadmesh(rtopo)
meshpts <- rgdal::project(t(qtopo$vb[1:2, ]), "+proj=lcc +lat_1=40 +lat_2=50 +lon_0=-60")
plot(raster::extent(meshpts), asp = 1)
for (i in seq_len(ncol(qtopo$ib))) {
  polygon(meshpts[qtopo$ib[,i], ], col = cols[i], border = NA)
}