...built with SveltR.

IxD:
Web:
Viz:

Real~Currents




 

Events & Shade Across Time & Terrain: Part 1

Given two texture images that describe the same terrain, one color and one grayscale, the first objective is to simulate the perceived time-of-day by altering the color texture with light and shade effects. The ray-tracing-based shade algorithm that I will use comes from the blog post Throwing Shade at the Cartographer Illuminati by Tyler Morgan-Wall.

shadow <- matrix(1, nrow = nrow(heightmap), ncol = ncol(heightmap))
tanangheight <- matrix(1, nrow = nrow(heightmap), ncol = ncol(heightmap))

sunangle <- 45 / 180 * pi
angle <- -90 / 180 * pi
diffangle <- 90 / 180 * pi
numberangles <- 25
anglebreaks <- seq(angle, diffangle, length.out = numberangles)
maxdistance <- floor(sqrt(nrow(heightmap)^2 + ncol(heightmap)^2))

for (y in 1:nrow(heightmap)) {
    for (x in 1:ncol(heightmap)) {
        for (anglei in anglebreaks) {
            for (k in 1:maxdistance) {
                xcoord <- (x + cos(sunangle) * k)
                ycoord <- (y + sin(sunangle) * k)
                if (xcoord > nrow(heightmap) ||
                    ycoord > ncol(heightmap) ||
                    xcoord < 0 ||
                    ycoord < 0) {
                    break
                } else {
                    tanangheight[y, x] <- heightmap[y, x] + tan(anglei) * k
                    if (all(
                        c(heightmap[ceiling(ycoord), ceiling(xcoord)],
                          heightmap[ceiling(ycoord), floor(xcoord)],
                          heightmap[floor(ycoord), ceiling(xcoord)],
                          heightmap[floor(ycoord), floor(xcoord)]
                        ) < tanangheight[y, x]
                    )) next
                    if (tanangheight[y, x] < bilinear(1:ncol(heightmap), 1:nrow(heightmap), heightmap, x0 = xcoord, y0 = ycoord)$z) {
                        shadow[y, x] <- shadow[y, x] - 1 / length(anglebreaks)
                        break
                    }
                }
            }
        }
    }
}

The R implementation of this particular method is pretty slow, as explained in the post about the algorithm referenced above, even on the small 128x128 example that was used. So, the author created a package called rayshader that imports a C implementation of the same. Here’s how it works:

elmat <- matrix(1, nrow = nrow(image_test), ncol = ncol(image_test))

for (y in 1:nrow(image_test)) {
    for (x in 1:ncol(image_test)) {
        elmat[x, y] <- heightmap[y, x]
    }
}

sunangle <- 45

elmat %>%
    ray_shade(
        maxsearch = floor(sqrt(nrow(elmat)^2 + ncol(elmat)^2)) / 64,
        zscale = 100, sunaltitude = 45, sunangle = sunangle, lambert = TRUE,
        multicore = TRUE) %>%
    plot_map()

Following the examples is pretty straight forward and incredibly useful once you understand the parameters. However, the actual data I want to render and superimpose on the colormap is quite large, so I will break it down into quadrants before processing:

xparts <- 2
yparts <- 2
width <- ncol(heights)
height <- nrow(heights)
quad_width <- floor(width / xparts)
quad_height <- floor(height / yparts)

#print(paste0("data width: ", width))
#print(paste0("data height: ", height))
#print(paste0("quad width: ", quad_width))
#print(paste0("quad height: ", quad_height))

image_shaded <- array(dim = c(height, width, 3))

#print(paste0("image width: ", ncol(image_shaded)))
#print(paste0("image height: ", nrow(image_shaded)))

for (j in 1:yparts) {
    for (i in 1:xparts) {
        quad_width_start <- ((i - 1) * quad_width)
        quad_width_end <- (i * quad_width)
        quad_height_start <- ((j - 1) * quad_height)
        quad_height_end <- (j * quad_height)
        elmat <- matrix(1, nrow = quad_height, ncol = quad_width)

        for (y in 1:quad_height) {
            for (x in 1:quad_width ) {
                elmat[x, y] <- heights[quad_height_start + y, quad_width_start + x] * height_mult
            }
        }

        sunangle <- 45

        revx_elmat <- matrix(1, nrow = nrow(elmat), ncol = ncol(elmat))
        for (y in 1:nrow(elmat)) {
            for (x in 1:ncol(elmat)) {
                revx_elmat[y, x] <- elmat[(nrow(elmat) - y) + 1, x]
            }
        }

        single_angle_shadow <- (revx_elmat %>%
            ray_shade(
                maxsearch = floor(sqrt(nrow(elmat)^2 + ncol(elmat)^2)) / 64,
                zscale = 100, sunaltitude = 45, sunangle = sunangle, lambert = TRUE,
                multicore = TRUE)
        )

        if (file.exists("data")) {
            rayshader_test_file <- paste0("data/terrain-rayshaded-", (i - 1), "", (j - 1), ".png")
        } else {
            rayshader_test_file <- file.choose(new = TRUE)
        }

        # image function and ggplot are slow; save with PNG function from PNG package.
        png::writePNG(t(single_angle_shadow), rayshader_test_file)

        cat(paste0('<img id="shademap', (i - 1), "", (j - 1), '" src="/post/', rayshader_test_file, '" style="width: 100%" />'))

        for (y in 1:quad_height) {
            for (x in 1:quad_width ) {
                # sanity check
                #image_shaded[quad_height_start + y, quad_width_start + x, ] <- c(
                #    image_color[quad_height_start + y, quad_width_start + x + (0 * width)],
                #    image_color[quad_height_start + y, quad_width_start + x + (1 * width)],
                #    image_color[quad_height_start + y, quad_width_start + x + (2 * width)]
                #)

                # combine shadow with colormap, with a value ratio of 1:2 respectively
                image_shaded[quad_height_start + y, quad_width_start + x, ] <- c(
                    ((single_angle_shadow[x, y] / 2 + 0.5) * image_color[quad_height_start + y, quad_width_start + x + (0 * width)]),
                    ((single_angle_shadow[x, y] / 2 + 0.5) * image_color[quad_height_start + y, quad_width_start + x + (1 * width)]),
                    ((single_angle_shadow[x, y] / 2 + 0.5) * image_color[quad_height_start + y, quad_width_start + x + (2 * width)])
                )
            }
        }
    }
}

if (file.exists("data")) {
    rayshaded_file <- paste0("data/terrain-shaded-colormap.png")
} else {
    rayshaded_file <- file.choose(new = TRUE)
}

## test write to png
png::writePNG(image_shaded, rayshaded_file)

cat(paste0('<img id="shademap" src="/post/', rayshaded_file, '" style="width: 100%" />'))

I have made use of the mapTerrainCoords method that was developed in Generating Planar Terrain Mesh to generate the terrain heights and normals, and then write them out to a PNG and JSON file, respectively: