Loading 3DEP elevation data in C++ with GDAL

Readme.md

Download ArcGrid 1/3 arc-second data from

https://viewer.nationalmap.gov/basic/?basemap=b1&category=ned&q=&zoom=10&bbox=-107.55203247,38.85788953,-106.02355957,39.64482504&preview=&avail=DEM 1/3 arc-second&refpoly=

(hit "find products", use "footprint" link to view on the map, then click "download")

We'll start with https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/13/ArcGrid/n40w107.zip.

Extract the zip file.

---

Install GDAL:

sudo apt install gdal-bin libgdal-dev

---

To do a quick sanity check of the data, gdalinfo will show you the bounding box, min and max elevations, etc.

$ gdalinfo grdn40w107_13 -stats

Driver: AIG/Arc/Info Binary Grid
Files: grdn40w107_13
       grdn40w107_13.aux.xml
       grdn40w107_13/w001001x.adf
       grdn40w107_13/prj.adf
       grdn40w107_13/dblbnd.adf
       grdn40w107_13/log
       grdn40w107_13/hdr.adf
       grdn40w107_13/sta.adf
       grdn40w107_13/w001001.adf
       grdn40w107_13/metadata.xml
Size is 10812, 10812
Coordinate System is:
GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9108"]],
    AUTHORITY["EPSG","4269"]]
Origin = (-107.000555555600002,40.000555555555515)
Pixel Size = (0.000092592592593,-0.000092592592593)
Corner Coordinates:
Upper Left  (-107.0005556,  40.0005556) (107d 0' 2.00"W, 40d 0' 2.00"N)
Lower Left  (-107.0005556,  38.9994444) (107d 0' 2.00"W, 38d59'58.00"N)
Upper Right (-105.9994444,  40.0005556) (105d59'58.00"W, 40d 0' 2.00"N)
Lower Right (-105.9994444,  38.9994444) (105d59'58.00"W, 38d59'58.00"N)
Center      (-106.5000000,  39.5000000) (106d30' 0.00"W, 39d30' 0.00"N)
Band 1 Block=256x4 Type=Float32, ColorInterp=Undefined
  Min=1885.109 Max=4397.064 
  Minimum=1885.109, Maximum=4397.064, Mean=3008.319, StdDev=471.554
  NoData Value=-3.4028234663852886e+38
  Metadata:
    STATISTICS_MAXIMUM=4397.064453125
    STATISTICS_MEAN=3008.318902968
    STATISTICS_MINIMUM=1885.109375
    STATISTICS_STDDEV=471.5536144229

---

Compile and run the code:

$ g++ -std=c++11 -Wall -Werror -g -O2 read_dem.cpp -lgdal -o read_dem
$ ./read_dem

Image is 10812 x 10812 px
Origin: -107.001, 40.0006 degrees
Pixel size: 9.25926e-05, -9.25926e-05 degrees
Loading array...done
39.1625,-106.97 maps to pixel 334.749,9051.18
Elevation there is 3771.7m

See http://www.gdal.org/gdal_tutorial.html for more info on the GDAL API

read_dem.cpp

#include <iostream>
#include <vector>
#include "gdal/gdal_priv.h"
#include "gdal/cpl_conv.h"

int main() {
    // Initialize GDAL
    GDALAllRegister();
    
    // Open the file
    GDALDataset  *dataset = (GDALDataset *) GDALOpen("grdn40w107_13", GA_ReadOnly);
    if( dataset == NULL ) {
      std::cout << "Failed to open" << std::endl;
      exit(1);
    }
    
    // Get image metadata
    unsigned width = dataset->GetRasterXSize();
    unsigned height = dataset->GetRasterYSize();
    
    std::cout << "Image is " << width << " x " << height << " px" << std::endl;
    
    double originX, originY, pixelSizeX, pixelSizeY;
    
    double geoTransform[6];
    if (dataset->GetGeoTransform(geoTransform) == CE_None ) {
        originX = geoTransform[0];
        originY = geoTransform[3];
        pixelSizeX = geoTransform[1];
        pixelSizeY = geoTransform[5];
    } else {
      std::cout << "Failed read geotransform" << std::endl;
      exit(1);
    }
    
    std::cout << "Origin: " << originX << ", " << originY << " degrees" << std::endl;
    std::cout << "Pixel size: " << pixelSizeX << ", " << pixelSizeY << " degrees" << std::endl; 
    
    // pixelSizeY is negative because the origin of the image is the north-east corner and positive
    // Y pixel coordinates go towards the south
    
    // Get the raster band (DEM has one raster band representing elevation)
    GDALRasterBand  *elevationBand = dataset->GetRasterBand(1);
    
    // Create an array of width*height 32-bit floats (~400MB memory)
    std::vector<float> data(width * height, 0.0f);
   
    // Read the entire file into the array (you can change the options to read only a portion
    // of the file, or even scale it down if you want)
     
    std::cout << "Loading array..." << std::flush;
    elevationBand->RasterIO(GF_Read, 0, 0, width, height, &data[0], width, height, GDT_Float32, 0, 0);
    std::cout << "done" << std::endl;
    
    // Close the file
    GDALClose(dataset); 
    
    // Be careful with axis order: coordinates are traditionally written in lat, lon order, but
    // those are Y, X, respectively.    
    double pointLat = 39.1624833;
    double pointLon = -106.9695603;
    
    double pixelX = (pointLon - originX) / pixelSizeX;
    double pixelY = (pointLat - originY) / pixelSizeY;
    
    std::cout << pointLat << "," << pointLon << " maps to pixel " << pixelX << "," << pixelY << std::endl;
    
    if (pixelX < 0 || pixelX > width || pixelY < 0 || pixelY > height) {
      std::cout << "Point out of bounds!" << std::endl;
      exit(1);
    }
    
    double elev = data[int(pixelX) + int(pixelY) * width];
    
    std::cout << "Elevation there is " << elev << "m" << std::endl;
}