Our work is done under the UN Vector Tile Toolkit activities and the UN Open GIS Initiative.
Introduction
When we develop a web map application using MapLibre GL JS or Mapbox GL JS, we often use elevation information for drawing hill shade and/or enabling 3D terrain expression. In such a case, elevation data is often treated as raster tile, specifically as RGB elevation tile that stores elevation value in RGB values. Therefore, I think that the RGB elevation tile is a kind of defact format for the elevation tile for web mapping.
On the other hand, however, it would be hard to find free and open RGB elevation tiles. If we want to use Mapbox Terrain-DEM or maptiler Terrain RGB from maptiler DATA, we have to pay for these data consumption. They host tile services and these are really useful for many users, but you can also choose hosting tiles by yourself by converting and hosting your elevation tile by yourself withouth any fee for data cousumption. In the world, there are several digital elevation model (DEM) available as open data such as SRTM by NASA/USGS and they can be used as the source of your RGB elevation tile.
As a part of our UN Vector Tile Toolkit activities under the UN Open GIS Initiative, I wanted to create RGB elevation tile by ourselves from open data (NASADEM) for our own purpose. This is the record of my work, and I hope this article may help those who would like to develop his/her own RGB elevation tile for free of charge.
(In August 2022, I created RGB tiles from SRTM, but it has a lot of void value due to the original data. This time, I will use NASADEM so that we have less void area.)
What is RGB elevation tile?
RGB elevation tile is raster tile that stores elevation infocation stored in RGB values. It is used for many mapping library such as Mapbox GL JS and MapLibre GL JS. We can create a 3D terrain map and/or hill shade map from it.
height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)
https://documentation.maptiler.com/hc/en-us/articles/4405444055313-RGB-Terrain-by-MapTiler
https://docs.mapbox.com/data/tilesets/reference/mapbox-terrain-rgb-v1/
Source data (NASADEM)
The Land Processes Distributed Active Archive Center (LP DAAC) is distributing NASADEM for the global coverage (note: no data in high lattitude area). Like other NASA/USGS's product, they are free of charge. Data can be downloaded from this URL: https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11/ (user account is needed.)
Figure 1. Comparison of data in Western Alps (n45-46, e007-008). Left: SRTM, Right: NASADEM
NASADEM is based on SRTM and comlemented by ASTER G-DEM and ALOS 30 meter DEM. It has 1 arc second spatial resolutions (about 30 meters), and has less void area compared with SRTM (Fig 1). A file is for 1 degree by 1 degree area with 3601 by 3601 pixels. There are more than 14 thousand files in total. For detail, we can refer to the following website.
My Environment
- Windows 10 Enterprise
- Docker v20.10.8
- PowerShell
- Google Chrome Version 107.0.5304.87 (64-bit)
Basic tactics
When we create the RGB elevation tile, converting for the range of source file (1 degree by 1 degree) each time would be tough because we needed to do 14,520 conversions and meeded to merge outputs for raster tiling. It would be also tough to convert the whole globe at once given the size of input file (we would need to merge all the source first...).
So, my stragegy or tactics, which is often used in UN Vector Tile Toolkit, is to use the spatial modules of the ZL6 (zoom level 6) tiles. We have 4096 tiles for ZL6, but about 2/3 of them are water area without DEM. In addition, the northern and southern limits of SRTM (and NASADEM) are lower than the web mercator projection limits. Therefore, we may need to do the work for less than thousand modules/spatial extents (i.e., it is 931 for NASADEM). The figure 6 is an example of ZL6 tiling over SRTM data range (green boxes, 1 by 1 degree, are SRTM bounding boxes). Please note this ZL6 tile extent is in the web meracator project, so some of them do not look like squre in the DD projection.
Figure 6. ZL6 tile extents are overlaied to SRTM data ranges (grenn cells, 1 degree by 1 degree, are SRTM data extent.)
Our basic processing workflow would be like this:
- For any (6(=z), x, y),
- Check the bounding box of that tile extent in decimal degree.
- Check if there is any NASADEM files in that extent. --> if not file, we do not need to work for that area.
- If there is/are NASADEM file(s), merger them as a single input.tif for the extent.
- Conver merged tiff file into RGB elevation tile in mbtiles format (a kind of SQlite format containing all tiles for that extent.)
Procedure
Step 1. Start Docker
In order to establish working environment, I ran a docker container. There is a docker image named unvt/rgbify-node for RGB elevation tile processing with nodejs.
docker run -it --rm -v ${PWD}:/data unvt/rgbify-node
cd /data
Step 2. Data download
As the LP DAAC of NASA provides data download via command line (curl), I downloded 14,520 files with the coomand line in about 20 hours. The following webpage helped me to get the data with curl command. I created my user ID for NASA EarthDATA at https://urs.earthdata.nasa.gov/users/new to get the data.
2-1 Create .netrc file
Because I used my user ID and password for NASA EarthDATA when we downloaded DEM via curl command, I created .netrc file first.
touch ~/.netrc | chmod og-rw ~/.netrc | echo machine urs.earthdata.nasa.gov >> ~/.netrc
echo login MYUSERNAME >> ~/.netrc
echo password MYPASSWORD >> ~/.netrc
You can refer to my another article (in Japanese) here for details: https://qiita.com/T-ubu/items/8f955f21a5603a3483b3
2-2 Data download using sh file
I prepared a shell file to download all 14520 files.
https://raw.githubusercontent.com/ubukawa/nasadem-rgbify/main/src/dl.sh
Figure 2. A screenshot of the sh file for download
Then, I ran the command to download the data.
mkdir src
cd src
./dl.sh
I downloaded all 14,520 files. Their data size in total is about 101GB.
Step 3. Unpack data
Because gdal_merge did not work with zipped hgt files, I unpacked them. All the zipped files were in the src directory, so I unpacked them with the following command. To save the storage, I deleted num file and swb file.
mkdir 01_unpack
for f in src/*.zip; do unzip ${f} -d 01_unpack; rm 01_unpack/*.num; rm 01_unpack/*.swb; done
Because there were 14,520 zipped files, I needed to wait for a while before unzip work started, and it took about two days to unpack all the files in my environment. I stored the data in the portable HDD, so it would be faster if I work with SSD storage.
Figure 3. A screenshot during the download
The total size of the unpacked file was more than 350GB. Each DEM is stored in 16 bit siged integer format with the extention of "hgt" meaning height.
GDAL and QGIS has a HGT driver, originally for SRTM HGT format, and they can deal with this format. Data itself has no information of the projection, spatial resolution and extent, etc., but georeference is done based on the filename that includes the min coordinate of x/y in decimal degree (e.g. "n27e027.hgt").
SRTMHGT – SRTM HGT Format: https://gdal.org/drivers/raster/srtmhgt.html
Step 4. Merge (at each ZL6 spatial module)
I merged all source DEM to prepare an input file for each spatial extent/module.
To avoid a memory overflow, I separated "merge" process from "rgbifying" process. (I tried doing them all at once, but I failed.)
The following note is for those who want to create RGB elevation for whole globe. You can refer to my another note if you want to work with a small number of source files: https://qiita.com/T-ubu/items/09720f25d5278d85222c
4-1. Checking the gdal_merge.py with hgt files
First, I tested some command and confirmed that gdal_merge.py work with hgt files. It confirmed that georeference information was also reflected to the output. (The test2.tig in the following figure has only 2 source files, but has about 75MB. Locations are well recognized.)
Figure 4. Testing the gdal_merge.py from hgt sources to create merged tiff
Each source was properly located when I merged.
Figure 5. Checking the output of the gdal_merge.py with QGIS (check the coordinates and blank in the center)
4-2. Preparing the script and running it
I have prepared some scripts. All files are stored here: https://github.com/ubukawa/nasadem-rgbify. In order to conduct the merge works for all spatial modules, I prepared:
nodejs package file (click here to open)
{
"name": "data",
"version": "1.0.0",
"description": "Rgbifying NASA",
"main": "index4merge.js",
"scripts": {
"merge": "node index4merge.js",
"rgbify": "node index4rgb.js"
},
"repository": {
"type": "git",
"url": "git+https://github.com/ubukawa/nasadem-rgbify.git"
},
"keywords": [],
"author": "",
"license": "ISC",
"bugs": {
"url": "https://github.com/ubukawa/nasadem-rgbify/issues"
},
"homepage": "https://github.com/ubukawa/nasadem-rgbify#readme",
"dependencies": {
"@mapbox/tilebelt": "^1.0.2",
"better-queue": "^3.8.10",
"child_process": "^1.0.2",
"config": "^3.3.7",
"fs": "^0.0.1-security",
"hjson": "^3.2.2"
}
}
config file in hjson (click here to open)
{
srcDir: 01_unpack
gdalmergePath: gdal_merge.py
rasterioPath: rasterio
maxZ: 12
minZ: 6
concurrent: 3
maxRetries: 3
retryDelay: 5000
mergeDir: merge
mbtilesDir: mbtiles
}
The congig file (config/default.hjson) containes not only parameter for merge work, but also parameters for rgbifying.
nodejs script for merge (click here to open)
//modules
const config = require('config')
const fs = require('fs')
const {spawn} = require('child_process')
const tilebelt = require('@mapbox/tilebelt')
const Queue = require('better-queue')
//config parameters
const srcDir = config.get('srcDir')
const mergeDir = config.get('mergeDir')
const mbtilesDir = config.get('mbtilesDir')
const gdalmergePath = config.get('gdalmergePath')
let modulesObj = {} //object {key: [srcFile, ... ], ...}
let emptyModules = []
let keys = [] //Array of key such as "6-x-y"
let keyInProgress = []
let idle = true
let countModule = 0
const isIdle = () => {
return idle
}
const sleep = (wait) => {
return new Promise((resolve, reject) => {
setTimeout(()=> {resolve(), wait})
})
}
let fileList = fs.readdirSync(srcDir) //list from the src folder (01_unpack)
let nasademFiles = fileList.filter(r => r.indexOf('.hgt') !== -1) //only hgt file
//keys (6-x-y)
for (x = 0; x < 64; x ++){
for (y = 0; y < 64; y++) {
let key = `6-${x}-${y}`
keys.push(key)
}
}
for (const key of keys){
//for (const key of ['6-31-31','6-32-32']){
let [tilez, tilex, tiley] = key.split('-')
tilex = Number(tilex)
tiley = Number(tiley)
tilez = Number(tilez)
const bbox = tilebelt.tileToBBOX([tilex, tiley, tilez])
modulesObj[key] = []
for (x=Math.floor(bbox[0]); x < bbox[2]; x++ ){
m = x.toString(10) // 10 means decimal
if(x < 0) {
m = m.replace("-","")
if(m.length == 1){
m = `00${m}`
} else if (m.length == 2) {
m = `0${m}`
}
m = `W${m}`
} else {
if(m.length == 1){
m = `00${m}`
} else if (m.length == 2) {
m = `0${m}`
}
m = `E${m}`
} // Then, m has proper string
for (y = Math.floor(bbox[1]); y < bbox[3]; y++){
n = y.toString(10)
if(y<0){
n = n.replace("-","")
if(n.length == 1) {
n = `0${n}`
}
n = `S${n}`
} else {
if(n.length == 1) {
n = `0${n}`
}
n = `N${n}`
}
nm = `${n.toLowerCase()}${m.toLowerCase()}.hgt` //e.g. e27e027.hgt
if(nasademFiles.includes(nm)){
//console.log (`${nm}---> yes(${key})`)
modulesObj[key].push(`${srcDir}/${nm}`)
}
}
}
if (Object.keys(modulesObj[key]).length == 0) {
emptyModules.push(key)
}
if (modulesObj[key].length == 0){
delete modulesObj[key]
}
}
const queue = new Queue(async (t, cb) => {
const key = t.key
const tile = t.tile
const [z, x, y] = tile
const mergedPath = `${mergeDir}/${key}.tif`
const tmpPath = `${mbtilesDir}/part-${key}.mbtiles`
const dstPath = `${mbtilesDir}/${key}.mbtiles`
countModule ++
keyInProgress.push(key)
console.log(`[${keyInProgress}] in progress`)
console.log(`--- ${key} (${countModule}/${Object.keys(modulesObj).length}) starts: ${modulesObj[key].length} src file/files`) //list of src files
//console.log(`--- ${key} (${countModule}/${Object.keys(modulesObj).length}): ${modulesObj[key].length} (${modulesObj[key]})`) //list of src files
let gdalmergeArray = [
'-o', mergedPath
]
const mgStartTime = new Date()
gdalmergeArray = gdalmergeArray.concat(modulesObj[key])
if(fs.existsSync(mergedPath)){
console.log(`--- ${key}: file already exists (${mgStartTime.toISOString()})`)
keyInProgress = keyInProgress.filter((v) => !(v === key))
return cb()
} else {
const gdalmerge = spawn(gdalmergePath, gdalmergeArray)
//gdalmerge.stdout.on('data', (data) => {
// console.log(`stdout: ${data}`)
//})
gdalmerge.stderr.on('data', (data) =>{
console.log(`stderr(at merge):${data}`)
})
gdalmerge.on('error', (error) => console.log(error.message.message))
gdalmerge.on('exit', (code, signal) =>{
if(code) console.log(`process exit with code: ${code}.`)
if(signal) console.log(`process killed with signal: ${signal}.`)
const mgEndTime = new Date()
console.log(`--- ${key}: ${modulesObj[key].length} src file/files merge ends (${mgStartTime.toISOString()} --> ${mgEndTime.toISOString()} )`)
keyInProgress = keyInProgress.filter((v) => !(v === key))
//fs.renameSync(tmpPath,dstPath)
//fs.unlinkSync(mergedPath)
return cb()
})
}
},{
concurrent: config.get('concurrent'),
maxRetries: config.get('maxRetries'),
retryDelay: config.get('retryDelay')
})
const queueTasks = () => {
for (let module of Object.keys(modulesObj)){
let tile = module.split('-').map(v => Number(v))
queue.push({
key: module,
tile: tile
})
}
}
const shutdown = () => {
console.log('** production system (merge) shutdown! (^_^) **')
}
const main = async () =>{
const stTime = new Date()
console.log(`-------UNVT---------------\n${stTime.toISOString()}: Production starts. \n- From the saved sources, we have ${Object.keys(modulesObj).length} modules with NASADEM. \n- ${emptyModules.length} modules are without NASADEM.\n- Here is the list of ${Object.keys(modulesObj).length} modules: \n${Object.keys(modulesObj)}\n--------------------------`)
queueTasks()
queue.on('drain', () => {
const closeTime = new Date()
console.log(`Production ends: ${stTime.toISOString()} --> ${closeTime.toISOString()}`)
shutdown()
})
}
main()
Then, I ran the following script.
npm install
npm run merge
It was turned out that 931 modules out of 4096 have NASADEM data in their extent. So, we will obtain 931 tiff files through this process. This statics was automatically calculated from the files in the directory where I stored the unpacked files.
Figure 6. A snapshot at the begining of the processing
For some modules, there are a few source files and it took short time. If a module has many files, it took longer time. Because I used the spatial modules based on the ZL6 tile, the maximum number of the source files in a module would be 49 and its output would be around a few GB.
Figure 7. A snapshot during merge work
It took about two days to merge files for all spatial modules that have NASADEM. The total size of all 931 merged files is 587GB. The total size increased compared with the total size of unpacked files, but it is natural because some source files (1 by 1 degrees) are duplicated for different ZL6 extents.
I remember the merge process took about 5 days when I worked with SRTM dem (source files are in tiff format) in August 2022. I changed the resource setting of the Docker since then, so it would affect the conversion time. Anyway, the log of the work can be accessible from here: https://github.com/ubukawa/nasadem-rgbify/blob/main/merge-log.txt
Figure 8. Merged files for respective spatial modules as inputs for RGBifying
Step 5. RGBifying
I prepared another nodejs script for RGBifying merged files.
nodejs script for rgbifying (click here to open)
//modules
const config = require('config')
const fs = require('fs')
const {spawn} = require('child_process')
const Queue = require('better-queue')
//config parameters
const mergeDir = config.get('mergeDir')
const mbtilesDir = config.get('mbtilesDir')
const rasterioPath = config.get('rasterioPath')
const maxZ = config.get('maxZ')
const minZ = config.get('minZ')
let keys = [] //Array of key such as "6-x-y"
let keyInProgress = []
let idle = true
let countModule = 0
const isIdle = () => {
return idle
}
const sleep = (wait) => {
return new Promise((resolve, reject) => {
setTimeout(()=> {resolve(), wait})
})
}
let mgfileList = fs.readdirSync(mergeDir) //list from the merge folder
mgfileList = mgfileList.filter(r => r.indexOf('.tif') !== -1) //only tiff file
for (let i=0; i<mgfileList.length; i++){
keys.push(mgfileList[i].replace('.tif',''))
}
const queue = new Queue(async (t, cb) => {
//const startTime = new Date()
const key = t.key
const tile = t.tile
const [z, x, y] = tile
const mergedPath = `${mergeDir}/${key}.tif`
const tmpPath = `${mbtilesDir}/part-${key}.mbtiles`
const dstPath = `${mbtilesDir}/${key}.mbtiles`
countModule ++
keyInProgress.push(key)
console.log(`[${keyInProgress}] in progress`)
console.log(`--- ${key} (${countModule}/${keys.length}) starts`) //list of src files
const rgbStartTime = new Date()
if(fs.existsSync(dstPath)){
console.log(`--- ${key}: file already exists (${rgbStartTime.toISOString()})`)
keyInProgress = keyInProgress.filter((v) => !(v === key))
return cb()
} else {
const rgbify = spawn(rasterioPath, [
'rgbify', '-b','-10000','-i','0.1', '--max-z', maxZ, '--min-z', minZ,
'--format', 'webp', '--bounding-tile', `[${x.toString()},${y.toString()},${z.toString()}]`,
mergedPath, tmpPath
])
//rgbify.stdout.on('data', (data) => {
// console.log(`stdout: ${data}`)
//})
rgbify.stderr.on('data', (data) =>{
console.log(`stderr(at rgbify):${data}`)
})
rgbify.on('error', (error) => console.log(error.message.message))
rgbify.on('exit', (code, signal) =>{
if(code) console.log(`process exit with code: ${code}.`)
if(signal) console.log(`process killed with signal: ${signal}.`)
keyInProgress = keyInProgress.filter((v) => !(v === key))
fs.renameSync(tmpPath,dstPath)
//fs.unlinkSync(mergedPath)
const rgbEndTime = new Date()
console.log(`--- ${key} ends: (${rgbStartTime.toISOString()} --> ${rgbEndTime.toISOString()} )`)
return cb()
})
}
},{
concurrent: config.get('concurrent'),
maxRetries: config.get('maxRetries'),
retryDelay: config.get('retryDelay')
})
const queueTasks = () => {
for (let key of keys){
//for (let tile of [[6,32,20],[6,32,21],[6,32,22],[6,32,23],[6,33,20],[6,33,21],[6,33,22]]){
//for (let key of ['bndl1', 'bndl2', 'bndl3', 'bndl4', 'bndl5', 'bndl6']){
let tile = key.split('-').map(v => Number(v))
queue.push({
key: key,
tile: tile
})
}
}
const shutdown = () => {
console.log('** production system shutdown! (^_^) **')
}
const main = async () =>{
const stTime = new Date()
console.log(`-------UNVT---------------\n${stTime.toISOString()}: Production starts. \n- From the merged files, we have ${keys.length} modules with SRTM DEM. \n- Here is the list of ${keys.length} modules: \n${keys}\n--------------------------`)
queueTasks()
queue.on('drain', () => {
const closeTime = new Date()
console.log(`Production ends: ${stTime.toISOString()} --> ${closeTime.toISOString()}`)
shutdown()
})
}
main()
Then, I ran the command as below:
npm run rgbify
Figure 9. A screenshot during rgbifying
I got the RGB elevations tiles in mbtiles format. It took about 2 days (41h28m) to obtain 931 mbtiles files. The total size is 184GB.
It is noteworthy that I got some future warning message during the rgbifying for 6-63-20, 6-63-21, .., 6-63-42. They are at the western limit of the web map mercator, and I thought it could be somehow related to the tile extent while I confirmed that their input files did not exist over the E180/W180 line.
Figure 10. A screenshot of future warning for 6-63-y
The zoom level range of the outputs is from ZL6 to ZL12. We can adjust the zoom level range by eiditing the config setting in config/default.hjson. Thinking of the resolution of the source data and the total size of outputs, I used 12 as the maxzoom.
Step 6. Hosting
Although the data hosting is not the main scope of this article, I would like to breifly explain how I host these RGB elevation tiles for our use.
Case 1: Hosting with nodejs/express server our hosting (Our case)
To save the stroage, I use "mbtiles" format, but it does not simply work with the static hosting because raster tile is generally requested as png format. We have established a simple web server using nodejs/express, and created a spetial module for raster tile (png) delivery from mbtiles with a npm module, called mapbox/mbtiles( npm package / GitHub Repository ). And, because we have 931 mbtiles, not a single file, we have adjusted routing so that we can properly respond to the RGB elevation tile requests. For detaile, you can check here: https://github.com/unvt/coesite3/blob/main/routes/rgbElev.js (and also check app.js in the upper level if needed).
Case 2: Hosting as static contents
By converting tiles from mbtiles format into png format, we can simply host them as static contents at any webserver.
unvt/rgbify-node has a tool called mbutil by mapbox, we can simply obtain png files by running the following command.
mb-util mbtiles/6-hoge-hoge.mbtiles zxy/6-hoge-hoge
From my experiece, I think total data size of png files is greater than that of mbtiles. Sometimes, it can be 3 or 4 times larger. So, if your data is over 100GB in mbtiles, it would be over 300GB in png. This is the reason why I chose to use mbtiles for my hosting.
Summary
Thus, I have explained how I converted RGB elevation tiles from NASADEM. NASADEM has less void compared with SRTM.
References