Skip to content
Permalink
Browse files

\added block diagram from DEM R script plus readme

  • Loading branch information...
brownag committed Mar 12, 2019
1 parent 68a9b37 commit d43d7755677be33cc70dcd8546379fe1c21ecead
@@ -0,0 +1,139 @@
# dem-to-block_diagram.R
# 03/11/2019
# @author: andrew brown;
# based on demo by dylan beaudette

library(rayshader)
library(rgl)
library(raster)
library(rgdal)
library(FedData)
library(viridis)
library(imager)

######
### SETUP
######

## 1. read shapefile for overlay (must cover full extent of elevation .TIF)
# for example, ssurgo data symbolized on musym
thematic_shp <- readOGR(dsn='L:/NRCS/MLRAShared/CA630/Archived_OFFICIAL_DB/final/FG_CA630_GIS_2018.gdb',
layer="ca630_a")

## 2. thematic attribute - the column name in shapefile attribute table
mu.col <- "MUSYM"

## 3. digital elevation model (TIFF, or other raster-compatible format) for a chunk of space
# e.g. pan to desired area in ArcMap, and Data > Export Data > By Data Frame
elev <- raster('lidar_Tm_test_5m.tif')
#elev <- raster('elev.tif')

## 4. OPTIONAL: resample raster input
## if needed, you can resample to some other resolution/grid size
## note that the DEM/hillshade and any derived overlays/shadows will be in this resolution

## copy elevation raster
# elev_template <- elev

## change raster resolution in template (leaving all else the same)
# res(elev_template) <- c(5, 5) #5 m by 5 m

## resample elevation raster to desired template
# elev <- resample(elev, elev_template)

## save to file
# writeRaster(elev, filename='lidar_Tm_test_5m.tif')

######
### END SETUP
######

# convert elevation raster -> matrix
elmat <- matrix(extract(elev, extent(elev), buffer=1000), nrow=ncol(elev), ncol=nrow(elev))

# calculate (rectangular) boundary of DEM, use that to cut the overlay shapefile

# note: there is no specific reason your extent polygon has to be rectangular
# but it is done here because rasters are commonly rectangular and we wan
extent.poly <- FedData::polygon_from_extent(elev)
extent.poly <- spTransform(extent.poly, CRS(proj4string(thematic_shp)))
thematic_shp <- crop(thematic_shp, y = extent.poly)

# assign numeric value that is 1:1 with mukey
thematic_shp$munum <- match(thematic_shp[[mu.col]], unique(thematic_shp[[mu.col]]))

# produce raster
thematic_shp <- spTransform(thematic_shp, CRS(proj4string(elev)))
theme <- rasterize(x = thematic_shp, y = elev, 'munum')

# inspect raster representation of theme musym
plot(theme)

# create RGB array from rasterized theme
tf <- "temp.png"
sgdf <- as(theme, 'SpatialGridDataFrame')
scl <- function(x) (x - min(na.omit(x))) / diff(range(na.omit(x)))

# OPTIONAL: set SPECIFIC COLORS for certain theme levels
#n.grp <- length(unique(values(theme)))
#first.colors <- viridis(n.grp)
#first.colors[5] <- "#000000" # set 5th symbol to black
## ETC.
#alt.cols <- col2rgb(first.colors[values(theme)])
#cols <- alt.cols

# select some colors that ideally span the color ramps, get RGB
cols <- col2rgb(rev(viridis(256))[scl(values(theme)) * 255 + 1])
# comment this line out if you are setting colors manually with alt.colors


# get the RGB channels and put them in the SGDF
sgdf$red <- cols[1,]
sgdf$grn <- cols[2,]
sgdf$blu <- cols[3,]

# write to PNG then read it back with imager... to get the color array
# this is convenient but theoretically just rescaling the RGB channels in the SGDF
# would be sufficient (i.e. no file output required per se)
writeGDAL(sgdf[c("red", "grn", "blu")], tf, type="Byte", mvFlag=255, drivername="PNG")
load.array <- as.array(load.image(tf))

# take just the R, G and B -- and shift around the dimensions a bit for rayshader
my.array <- load.array[,,1,1:3]

# compute shadows
raymat <- ray_shade(elmat)
ambmat <- ambient_shade(elmat)

# with big files save the intermediates in case rgl crashes R or something
#save(raymat, ambmat, file = "intermediates.Rda")

# example: add overlay to static map
# elmat %>%
# sphere_shade(texture = "desert") %>%
# add_overlay(my.array, alphacolor=1) %>%
# plot_map()

# INTERACTIVE 3D PLOT WITH RGL

# set perspective with right-mouse + drag
# zoom with mouse wheel
# rotate with left-mouse + drag

# interactive 3D plot via rgl
elmat %>%
sphere_shade(texture = "desert") %>%
add_overlay(my.array, alphacolor=1) %>%
#add_water(detect_water(elmat), color="desert") %>%
add_shadow(raymat, max_darken = 0.4) %>%
add_shadow(ambmat, max_darken = 0.4) %>%
plot_3d(elmat, zscale=0.9, fov=0, theta=30, zoom=0.75, phi=45, windowsize = c(1000,800), lineantialias = TRUE)

# take a static picture of the rgl window
render_snapshot()

# important to clear the previous rgl window if any settings are adjusted
rgl::rgl.clear()

# # not well supported on Windows
# render_label(elmat, x=100, y=100, z=4000, zscale=20, text = "Somethign", textsize = 2, linewidth = 5, freetype=FALSE, antialias = TRUE)
@@ -0,0 +1,102 @@
---
title: "DEM to Block Diagram - Rayshader Demo"
author: "Andrew Brown"
date: "Last updated: March 12, 2019"
output:
html_document:
toc: true
toc_float: true
number_sections: true
---

# Get the script

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

This document explains how to use the `rayshader` R package to create 3-D thematic block diagrams from digital elevation models and a thematic shapefile.

Download the __R__ script here: [dem-to-block_diagram.R](dem-to-block_diagram.R)

# Install dependencies

If you don't already have the necessary packages, install them:
```{r, eval=F}
install.packages('rayshader','rgl','raster','rgdal','imager','FedData','viridis')
#note: FedData is only used for creating polygon from rectangular extent
#note: viridis is only used to generate a colorblind friendly color scheme
```

# Setup

There are a few options you will need to customize to run the script with your own data.

## Thematic Shapefile

Set the path to the shapefile/feature class you want to overlay on your 3D landscape.

```{r, eval=F}
thematic_shp <- readOGR(dsn = 'L:/PATH/TO/GEODATABASE/OR/FOLDER.gdb',
layer = "yourfeatureclass")
```

## Thematic Attribute
Set the _thematic attribute_ -- the column name in shapefile attribute table containing what you would like to symbolize.

```{r, eval=F}
mu.col <- "MUSYM"
```

## Digital Elevation Model
Load a raster digital elevation model (TIFF, or other raster-compatible format) for a chunk of space. You can easily create this for a desired extent by panning to area in ArcMap, and _Data_ > _Export Data_ > _By Data Frame_.

```{r, eval=F}
elev <- raster('C:/path/to/your_dem.tif')
```

### OPTIONAL: Resample Raster Input

If needed, you can resample your elevation raster to some other resolution/grid size. Note that the DEM/hillshade and any derived overlays/shadows will all be in the _same resolution_ as the elevation matrix.

If you have very detailed/large rasters, creating the elevation derivatives will take a long time and _rgl_ (3D visualization package) will run slow.

```{r, eval=F}
## copy elevation raster
elev_template <- elev
## change raster resolution in template (leaving all else the same)
res(elev_template) <- c(5, 5) #5 m by 5 m cells (from 1m x 1m)
## resample elevation raster to desired template
elev <- resample(elev, elev_template)
## save to file
writeRaster(elev, filename='lidar_resampled_5m.tif')
```

## Setting Custom Colors
After the input thematic shapefile has been rasterized, _viridis_ colors are assigned. By default, a set of colors that spans the color ramp that matches the number of unique levels in `mu.col` is used to render the map.

There is a section of commented-out code that allows you to modify the default set of viridis colors. Here is the code that you would need to edit.

```{r, eval=F}
# OPTIONAL: set SPECIFIC COLORS for certain theme levels
# calculate number of groups
n.grp <- length(unique(values(theme)))
# generate color palette with n.grp colors
first.colors <- viridis(n.grp)
# set 5th symbol to black (for example)
first.colors[5] <- "#000000"
## ETC.
alt.cols <- col2rgb(first.colors[values(theme)])
# remove/comment out other code to create `cols` (vector of colors)...
cols <- alt.cols
```

Large diffs are not rendered by default.

Binary file not shown.
BIN +81.9 KB BlockDiagram/temp.png
Binary file not shown.

0 comments on commit d43d775

Please sign in to comment.
You can’t perform that action at this time.