Previous tutorial: Adding a track bar to our application
Next tutorial: Video input using OpenCV and similarity measurement
| Original author | Marvin Smith |
|---|---|
| Compatibility | OpenCV >= 3.0 |
Geospatial raster data is a product heavily used in geographic information systems and photogrammetry. Raster data can often represent images and digital elevation models (DEMs). The standard library for loading GIS images is the Geographic Data Abstraction Library (GDAL). In this example, we will demonstrate techniques for loading a GIS raster format using OpenCV native functions. Additionally, we will provide examples of how OpenCV can use this data for new and interesting purposes.
Goals
The main goal of this tutorial is to
- How to load satellite images using OpenCV imread.
- How to use OpenCV imread – Load SRTM digital elevation model
- Given an image and the angular coordinates of a DEM, associate the elevation data with the image to find the elevation of each pixel.
- Shows a simple and easy to use terrain heat map example.
- Demonstrates the basic use of DEM data with orthophotos.
To achieve these goals, the following code takes as input a digital elevation model and a GeoTiff image of San Francisco. After processing the image and DEM data, a topographic heat map of the image is generated and identifies areas of the city that would be affected when the bay water level rises by 10, 50 and 100 meters.
Code
/*
* gdal_image.cpp -- Load GIS data into an OpenCV container using the geospatial data abstraction library
*/
// OpenCV header file
#include "opencv2/core.hpp
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp" */ // C++ standard library
// C++ standard library
#include <cmath
#include <iostream
#include <stdexcept
#include <vector
using namespace std;
//define corner points
// Note that the GDAL library can determine this natively
cv::Point2d tl(-122.441017, 37.815664);
cv::Point2d tr(-122.370919, 37.815311);
cv::Point2d bl(-122.441533, 37.747167);
cv::Point2d br(-122.3715, 37.746814);
// Determine the presentation angle
cv::Point2d dem_bl( -122.0, 38);
cv::Point2d dem_tr( -123.0, 37);
//Heat map color range
std::vector<std::pair<cv::Vec3b,double> > color_range;
// List of all function prototypes
cv::Point2d lerp( const cv::Point2d & amp;, const cv::Point2d & amp;, const double & amp; );
cv::Vec3b get_dem_color( const double & amp; );
cv::Point2d world2dem( const cv::Point2d & amp;, const cv::Size & amp;);
cv::Point2d pixel2world( const int & amp;, const int & amp;, const cv::Size & amp; );
void add_color( cv::Vec3b & amp; pix, const uchar & amp; b, const uchar & amp; g, const uchar & amp; r );
/*
* Linear interpolation
* p1 - point 1
* p2 - point 2
*t - Ratio from point 1 to point 2
*/
cv::Point2d lerp( cv::Point2d const & amp; p1, cv::Point2d const & amp; p2, const double & amp; t ){<!-- -->
return cv::Point2d( ((1-t)*p1.x) + (t*p2.x)、
((1-t)*p1.y) + (t*p2.y));
}
/*
* Interpolation color
*/
template <typename DATATYPE, int N>
cv::Vec<DATATYPE,N> lerp( cv::Vec<DATATYPE,N> const & minColor、
cv::Vec<DATATYPE,N> const & maxColor、
double const & amp; t ){<!-- -->
cv::Vec<DATATYPE,N> output;
for( int i=0; i<N; i + + ){<!-- -->
output[i] = (uchar)(((1-t)*minColor[i]) + (t * maxColor[i]));
}
return output;
}
/*
* Calculate demo color
*/
cv::Vec3b get_dem_color( const double & amp; elevation ){<!-- -->
// If the altitude is lower than the minimum value, return the minimum value
if( elevation < color_range[0].second ){<!-- -->
return color_range[0].first;
}
// If the altitude is higher than the maximum value, return the maximum value
if( elevation > color_range.back().second ){<!-- -->
return color_range.back().first;
}
// Otherwise, find a suitable starting index
int idx=0;
double t = 0;
for( int x=0; x<(int)(color_range.size()-1); x + + ){<!-- -->
// If the current height is lower than the next item, use the current
// Two colors as our range
if( elevation < color_range[x + 1].second ){<!-- -->
idx=x;
t = (color_range[x + 1].second - elevation)/ (color_range[x + 1].second - elevation)/ (color_range[x + 1].second)
(color_range[x + 1].second - color_range[x].second);
break;
}
}
// interpolate color
return lerp(color_range[idx].first, color_range[idx + 1].first, t);
}
/*
* Given pixel coordinates and the size of the input image, calculate the position of the pixel on the DEM image.
* Pixel location on the DEM image.
*/
cv::Point2d world2dem( cv::Point2d const & amp; coordinate, const cv::Size & amp; dem_size ){<!-- --> // Associate it with the DEM point.
// Associate it with the DEM point
// Assume that dem data is orthogonally corrected
double demRatioX = ((dem_tr.x - coordinate.x)/(dem_tr.x - dem_bl.x));
double demRatioY = 1-((dem_tr.y-coordinate.y)/(dem_tr.y-dem_bl.y));
cv::Point2d output;
output.x = demRatioX * dem_size.width;
output.y = demRatioY * dem_size.height;
return output;
}
/*
* Convert pixel coordinates to world coordinates
*/
cv::Point2d pixel2world( const int & amp; x, const int & amp; y, const cv::Size & amp; size ) {<!-- --> // Calculate the ratio of pixel position to size.
// Calculate the ratio of pixel position to its size
double rx = (double)x / size.width;
double ry = (double)y / size.height;
// Calculate LERP for each coordinate
cv::Point2d rightSide = lerp(tr, br, ry);
cv::Point2d leftSide = lerp(tl, bl, ry);
// Calculate the actual latitude/longitude coordinates of the interpolated coordinates
return lerp(leftSide, rightSide, rx);
}
/*
* Add color to specific pixel color value
*/
void add_color( cv::Vec3b & amp; pix, const uchar & amp; b, const uchar & amp; g, const uchar & amp; r ){<!-- -->
if( pix[0] + b < 255 & amp; & amp; pix[0] + b >= 0 ) {<!-- --> pix[0] + = b; }
if( pix[1] + g < 255 & amp; & amp; pix[1] + g >= 0 ){<!-- --> pix[1] + = g; }
if( pix[2] + r < 255 & amp; & amp; pix[2] + r >= 0 ){<!-- --> pix[2] + = r; }
}
/*
* Main function
*/
int main( int argc, char* argv[] ){<!-- -->
/*
* Check input parameters
*/
if( argc < 3 ){<!-- -->
cout << "usage: " << argv[0] << " <image_name> <dem_model_name>" << endl;
return -1;
}
// Load the image (note that we have no projection information. You will
// Need to load it yourself or use the full GDAL driver. These values are predefined
// at the top of this file
cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR);
//Load dem model
cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
//Create output product
cv::Mat output_dem( image.size(), CV_8UC3 );
cv::Mat output_dem_flood( image.size(), CV_8UC3 );
// For sanity's sake, make sure GDAL loads this as a signed short
if( dem.type() != CV_16SC1 ){<!-- --> throw std::runtime_error("DEM image type must be CV_16SC1"); } }
// Define the color range to create the output DEM heatmap
// Pair format ( Color, elevation ); Push from low to high
// NOTE: This is great for config files, but here is for a working demo.
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 188, 154, 46), -1));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 110, 220, 110), 0.25));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 150, 250, 230), 20);
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 160, 220, 200), 75));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 220, 190, 170), 100));
color_range.push_back( std::pair<cv::Vec3b,double>(cv::Vec3b( 250, 180, 140), 200));
//Define the minimum elevation
double minElevation = -10;
// Traverse each pixel in the image and calculate the presentation point
for( int y=0; y<image.rows; y + + ){<!-- -->
for( int x=0; x<image.cols; x + + ){<!-- -->
//Convert pixel coordinates to latitude/latitude coordinates
cv::Point2d coordinate = pixel2world(x, y, image.size());
// Calculate dem image pixel coordinates based on latitude/longitudinal coordinates
cv::Point2d dem_coordinate = world2dem( coordinate, dem.size() );
//Extract altitude
doubledz;
if( dem_coordinate.x >= 0 & amp; & amp; dem_coordinate.y >= 0 & amp; & amp;
dem_coordinate.x < dem.cols & amp; & amp; dem_coordinate.y < dem.rows ){<!-- -->
dz = dem.at<short>(dem_coordinate);
}else{<!-- -->
dz = minElevation;
}
//Write pixel values to file
output_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x);
// Calculate the color of the heatmap output
cv::Vec3b actualColor = get_dem_color(dz);
output_dem.at<cv::Vec3b>(y,x) = actualColor;
// Show the impact of 10 meters of sea level rise
if( dz < 10 ){<!-- -->
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 90, 0, 0);
}
// Show the impact of 50 meters of sea level rise
else if( dz < 50 ){<!-- -->
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 90, 0);
}
// Show the impact of 100 meters of sea level rise
else if( dz < 100 ){<!-- -->
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 0, 90);
}
}}
// print heat map
cv::imwrite( "heat-map.jpg" , output_dem );
//Print water immersion effect image
cv::imwrite( "flooded.jpg", output_dem_flood);
return 0;
}
How to read raster data using GDAL
This demonstration uses the default OpenCV imread function. The main difference is that in order to force GDAL to load an image, the appropriate flags must be used.
cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR);
When loading a digital elevation model, the actual value of each pixel is critical and cannot be scaled or truncated. For example, in image data, a pixel represented as a binary number with a value of 1 has the same appearance as a pixel represented as an unsigned character with a value of 255. For terrain data, the pixel value represents the altitude in meters. To ensure that OpenCV preserves the original values, use the GDAL flag and the ANYDEPTH flag in imread.
//Load dem model cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
If the type of the DEM model to be loaded is known in advance, it may be safer to test Mat::type() or Mat::depth() using assertions or other mechanisms. NASA or DOD specification documents can provide input types for various elevation models. The main types SRTM and DTED are signed dashes.
Comments
Latitude/longitude (geographical) coordinates should generally be avoided
Geographic coordinate systems are a spherical coordinate system, which means it is technically incorrect to use them in Cartesian mathematics. This demonstration uses these coordinate systems for readability and to be sufficiently accurate. A better coordinate system is the Universal Transverse Mercator coordinate system.
Find corner coordinates
An easy way to find the corner coordinates of an image is to use the command line tool gdalinfo. For orthorectified images that include projection information, use USGS EarthExplorer.
\f$> gdalinfo N37W123.hgt Driver: SRTMHGT/SRTMHGT File Format Files: N37W123.hgt Size is 3601, 3601 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", ... more output ... Corner Coordinates: Upper Left (-123.0001389, 38.0001389) (123d 0' 0.50"W, 38d 0' 0.50"N) Lower Left (-123.0001389, 36.9998611) (123d 0' 0.50"W, 36d59'59.50"N) Upper Right (-121.9998611, 38.0001389) (121d59'59.50"W, 38d 0' 0.50"N) Lower Right (-121.9998611, 36.9998611) (121d59'59.50"W, 36d59'59.50"N) Center (-122.5000000, 37.5000000) (122d30' 0.00"W, 37d30' 0.00"N) ... more output ...
Results
The following is the output of the program. Use the first image as input. For the DEM model, download the SRTM file located at USGS here. http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip
Input image
Heat Map
Heat map overlay