By Jay Bowen - GIS Specialist, The Digital Scholarship & Publishing Studio, The University of Iowa Libraries
The BTAA Geoportal has a wealth of interesting and beautiful historic scanned maps available to the public to download. Recently, I discovered John Henry Renshawe's shaded relief maps of US national parks from the early 1900s. Using his Panoramic View of the Yosemite National Park, California from 1914, I wanted to demonstrate an open-source technique for adding three dimensionality to these fantastic relief maps.
Download and Georeference the JPG File
The first step is to download the map here and load it into a georeferencer tool in your favorite GIS software. While I find georeferencing in ArcGIS Pro to be an intuitive breeze, I used the Georeferencer tool in QGIS to keep with an open-source and MacOS-oriented workflow. You can do a lot of amazing GIS work in the comfort of your home with a MacBook!
As shown in the screenshot above, I rubbersheeted the scanned Yosemite map by locating a series of six control points shared between the scanned map in the Georeferencer and the Esri Topo World basemap in the QGIS map window. If you are lucky, the map you wish to georeference will have latitude and longitude lines, and you can enter more accurate coordinates for your points at their intersections. However, in this case, I had to locate a few prominent shared features. For my transformation settings, I used a thin plate spline transformation type and a cubic spline resampling method. I find these settings best for older maps, especially when the projection is unknown.
Download Raster Elevation Data, Merge, and Clip
Next, you will need some raster elevation data at an
appropriate resolution for the area of interest. This will be used to generate
the three dimensional terrain used to warp the shaded relief map. I downloaded
two hgt files from a collection called "NASA
Shuttle Radar Topography Mission Global 1 arc second V003" at NASA's Earthdata
Search as shown in the image below.
Then, drag and drop the two raster hgt files into QGIS. You will need to stitch these two files together using the Merge tool, located at Raster > Miscellaneous > Merge within the dropdown options at the top. The merged image will exceed the extent of the shaded relief map, but you can clip the result by creating and using a mask layer. Choose Layer > Create Layer > New Shapefile Layer from the dropdown options at the top and trace the border of the mapped region, as shown below.
Now, Open Raster > Extraction > Clip Raster by Mask Layer. Enter the merged elevation raster as the input layer and the newly created shapefile as the mask layer. Run the tool. The result should look something like the following screenshot. You may notice a black edge in the output elevation raster overlapping the frame bordering the map. This is fine because these cells carry an elevation value of 0 and this part of the resulting three dimensional map will remain flat.
Define a Color and Elevation Range and Create an RGB.tif
You will need to convert the elevation raster to a colorized
tif file. This is not difficult, but it does require you to have gdal
installed. If you do not have this installed, you can use Homebrew
on a Mac to do so by pasting and executing the following command in
Terminal:
brew install gdal
Next, in Terminal, cd to the folder containing your clipped raster elevation file and create a txt file called color_relief.txt with the following command:
touch color_relief.txt
Now, open up this file and enter the following text:
3988 253 231 37
Finally, you need to return to Terminal to run the following gdaldem command on your clipped elevation raster. Make sure you cd to the path of the folder holding your raster and change "input_dem.tif" and "output_rgb.tif" to the names of your own files.
gdaldem color-relief input_dem.tif color_relief.txt output_rgb.tif
If you executed the command correctly, you should see an output file that looks like the image below.
Create PMTiles
Now that you have both a georeferenced map and an rgb tif,
you need to generate pmtiles from each of these for use with a MapLibre GL JS
web map application. First, before you generate the pmtiles, you should create
mbtiles. To do so, you will need rio-mbtiles.
Open Terminal in the folder holding your data. For the rgb tif, use the
following command to generate mbtiles:
rio mbtiles output_rgb.tif output_rgb.mbtiles --format PNG --zoom-levels 0..14 --tile-size 256 --resampling bilinear
Next, you will need the pmtiles command line tool. To convert the output mbtiles to pmtiles, just use the following command in Terminal:
pmtiles convert output_rgb.mbtiles output_rgb.pmtiles
You can repeat this process for the georeferenced relief map.
Code the Application with MapLibre GL JS
Now that you have all the files needed to generate the map,
you will want to set up your project folder with an "index.html" file
and a subdirectory for the pmtiles.
Open the index.html file and paste the code below. You can read the comments within the code to find out more about what each section is doing. If your file names are different from mine, make sure to change these in the code.
<html lang="en">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Yosemite National Park</title>
<script src="https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.js"></script>
<link href="https://unpkg.com/maplibre-gl@2.4.0/dist/maplibre-gl.css" rel="stylesheet" />
<script src="https://unpkg.com/pmtiles@2.5.0/dist/index.js"></script>
<style>
/* Define dimensions for the body of the page */
body {
margin: 0;
padding: 0;
}
#map {
position: absolute;
top: 0;
bottom: 0;
width: 100%;
}
</style>
</head>
<body>
<div id="map"></div>
<script>
let protocol = new pmtiles.Protocol();
const map = new maplibregl.Map({
style: {
version: 8,
sources: {
/* The topo source is the pmtiles for the scanned map downloaded from the geoportal. You can name these sources anything as long as your names are consistent throughout. */
topo: {
type: "raster",
url: "pmtiles://tiles/yosemite_topo.pmtiles", /* make sure to use your own file name here */
tileSize: 256,
minzoom: 0,
maxzoom: 14,
},
// The terrain source should be the rgb elevation tiles
terrainSource: {
type: "raster-dem",
url: "pmtiles://tiles/yosemite_rgb.pmtiles", /* make sure to use your own file name here */
tileSize: 256,
}
},
/* In this case, there is only one layer, which uses the topo source identified previously */
{
id: "topo",
type: "raster",
source: "topo",
}
],
/* The terrain uses the terrain source identified previously */
terrain: {
source: "terrainSource",
exaggeration: 0.005, /* adjust the exaggeration to your liking */
},
},
center: [-119.481452, 37.75339], /* add the center coordinates for your data */
zoom: 10, /* the initial zoom */
pitch: 40, /* the initial pitch */
bearing: 0, /* the initial bearing */
maxPitch: 85, /* the maximum allowed pitch */
maxZoom: 14 /* the maximum allowed zoom */
});
/* Add a navigation control to adjust zoom, pitch, and bearing */
map.addControl(
new maplibregl.NavigationControl({
showZoom: true,
showCompass: true,
})
);
/* Add a terrain control to allow users to toggle the 3D terrain */
map.addControl(
new maplibregl.TerrainControl({
exaggeration: 0.005, /* match to the exaggeration defined previously */
})
);
</script>
</body>
</html>
After this, you should have a working interactive 3D
terrain map of the scanned historic map of Yosemite National Park
downloaded from the BTAA
Geoportal using the MapLibre
GL JS TypeScript library. Enjoy your 3D mapping!
