Relief

Repo

Above Ground Shaded Reliefs

Join along @boyd.xyz/w

Yesterday in the Great Smokies
Webcam archive

Using lidar derivatives in ArcGIS Pro, Blender, and Python Notebooks

Boyd Shearer

Lecturer, UKy Department of Geography.
New Maps Plus, aka Digital Mapping.
outrageGIS mapping topographic trail maps.
Contact.

Long tradition in cartography of making shaded reliefs.

There's a ArcGIS Pro style for this.

Image sources

Erwin Raisz: Landforms of the United States, 1957
Tom Patterson: Equal Earth, 2018
Eduard Imhof: NE Switzerland, 1947, Mt. Everest, 1962
Swisstopo: 2021
Google Maps: 2022

Accuracy & Atmosphere

LiDAR -> above-ground

Hillshade a DSM

Blender 3D

Professional, free and open source 3D software
enables ray tracing to create subtle shading effects.

Ray tracing

Renders an atmospheric quality
by illuminating nearby objects with reflected light.
Takes a long minute to render.

How to use in a GIS?

Examples

Raster tile sets
Michele McHugh's geospatial PDF of Wolf Creek Dam

Methods

Workflow

Create DSM in ArcGIS Pro.
- Ky_DSM_First_Return_5FT_Phase1
Export as scaled 16-bit TIFF raster with world file.
Render in Blender matching original raster aspect ratio.
Using world file, bring raster back into AP and color with NDVI.

Python Notebook

Some notable cells from the notebook...

```py [1-2|4-5|7-8]
# Above ground, first return DSM for entire state
dsm = 'https://kyraster.ky.gov/arcgis/rest/services/ElevationServices/Ky_DSM_First_Return_5FT_Phase1/ImageServer'

# Load raster
a = arcpy.sa.Raster(dsm)

# verify
a
```
___
```py [1-3|5-6|8-11]
# Extract AOI and save to new layer
b = arcpy.sa.ExtractByMask(a, 'square_aoi')
b.save('original')

# Copy original raster to a 16-bit integer raster.
arcpy.management.CopyRaster(b, 'temp', '#',  '#',  '#',  '#',  '#',  '16_BIT_UNSIGNED', True)

# Load raster and get statistics
c = arcpy.sa.Raster('temp')
c.getStatistics()
```
___
```py [1-3|5-8]
# load low and high values
lo = c.minimum
hi = c.maximum

# Rescale raster to full range of integer values in 16-bit format, 0-65,535
d = ((c - lo) / (hi - lo))*65535
print(d.getStatistics())
d.save('scaled')
```
___
```py [1-2|4-6]
# Export TIFF
arcpy.management.CopyRaster('scaled', f'{root}\\{folder}\\elevations.tif', '#',  '#',  '#',  '#',  '#',  '16_BIT_UNSIGNED', True)

# Get pixel dimensions of raster for setting Blender parameters
size = arcpy.GetRasterProperties_management(d, 'COLUMNCOUNT')
print(f'Use this dimension for the X, Y size in Blender: {size.getOutput(0)}')
```
___
```py [1-4|1-7]
# Get aerial imagery
naip = 'https://kyraster.ky.gov/arcgis/rest/services/ImageryServices/Ky_NAIP_2020_2ft/ImageServer'
e = arcpy.sa.Raster(naip)
f = arcpy.sa.ExtractByMask(e, 'square_aoi')

#Verify
f
```
___
```py [1-3|1-6]
# Make and save NDVI
g = arcpy.sa.NDVI(f)
g.save('ndvi')

# Verify
g
```

Blender time

Make sure you have a 16-bit TIFF with a world file
and you know it's pixel dimensions.
Use this Blender project file
with the plane, camera, and light already set up.