Prev Tutorial: Adding a Trackbar to our applications!
Next Tutorial: Video Input with OpenCV and similarity measurement
| |
Original author | Marvin Smith |
Compatibility | OpenCV >= 3.0 |
Geospatial raster data is a heavily used product in Geographic Information Systems and Photogrammetry. Raster data typically can represent imagery and Digital Elevation Models (DEM). The standard library for loading GIS imagery is the Geographic Data Abstraction Library (GDAL). In this example, we will show techniques for loading GIS raster formats using native OpenCV functions. In addition, we will show some an example of how OpenCV can use this data for novel and interesting purposes.
Goals
The primary objectives for this tutorial:
- How to use OpenCV imread to load satellite imagery.
- How to use OpenCV imread to load SRTM Digital Elevation Models
- Given the corner coordinates of both the image and DEM, correlate the elevation data to the image to find elevations for each pixel.
- Show a basic, easy-to-implement example of a terrain heat map.
- Show a basic use of DEM data coupled with ortho-rectified imagery.
To implement these goals, the following code takes a Digital Elevation Model as well as a GeoTiff image of San Francisco as input. The image and DEM data is processed and generates a terrain heat map of the image as well as labels areas of the city which would be affected should the water level of the bay rise 10, 50, and 100 meters.
Code
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
std::vector<std::pair<cv::Vec3b,double> > color_range;
((1-t)*p1.
y) + (t*p2.
y));
}
template <typename DATATYPE, int N>
double const& t ){
for( int i=0; i<N; i++ ){
output[i] = (
uchar)(((1-t)*minColor[i]) + (t * maxColor[i]));
}
return output;
}
cv::Vec3b get_dem_color(
const double& elevation ){
if( elevation < color_range[0].second ){
return color_range[0].first;
}
if( elevation > color_range.back().second ){
return color_range.back().first;
}
int idx=0;
double t = 0;
for( int x=0; x<(int)(color_range.size()-1); x++ ){
if( elevation < color_range[x+1].second ){
idx=x;
t = (color_range[x+1].second - elevation)/
(color_range[x+1].second - color_range[x].second);
break;
}
}
return lerp( color_range[idx].first, color_range[idx+1].first, t);
}
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));
output.
x = demRatioX * dem_size.
width;
output.
y = demRatioY * dem_size.
height;
return output;
}
double rx = (double)x /
size.width;
double ry = (double)y /
size.height;
return lerp( leftSide, rightSide, rx );
}
if( pix[0] + b < 255 && pix[0] + b >= 0 ){ pix[0] += b; }
if( pix[1] + g < 255 && pix[1] + g >= 0 ){ pix[1] += g; }
if( pix[2] + r < 255 && pix[2] + r >= 0 ){ pix[2] += r; }
}
int main(
int argc,
char* argv[] ){
if( argc < 3 ){
cout << "usage: " << argv[0] << " <image_name> <dem_model_name>" << endl;
return -1;
}
if( dem.
type() !=
CV_16SC1 ){
throw std::runtime_error(
"DEM image type must be CV_16SC1"); }
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));
double minElevation = -10;
for(
int y=0; y<image.
rows; y++ ){
for(
int x=0; x<image.
cols; x++ ){
double dz;
if( dem_coordinate.
x >= 0 && dem_coordinate.
y >= 0 &&
dem_coordinate.
x < dem.
cols && dem_coordinate.
y < dem.
rows ){
dz = dem.
at<
short>(dem_coordinate);
}else{
dz = minElevation;
}
if( dz < 10 ){
add_color( output_dem_flood.at<
cv::Vec3b>(y,x), 90, 0, 0 );
}
else if( dz < 50 ){
add_color( output_dem_flood.at<
cv::Vec3b>(y,x), 0, 90, 0 );
}
else if( dz < 100 ){
add_color( output_dem_flood.at<
cv::Vec3b>(y,x), 0, 0, 90 );
}
}}
return 0;
}
n-dimensional dense array class
Definition mat.hpp:812
MatSize size
Definition mat.hpp:2160
_Tp & at(int i0=0)
Returns a reference to the specified array element.
int cols
Definition mat.hpp:2138
int rows
the number of rows and columns or (-1, -1) when the matrix has more than 2 dimensions
Definition mat.hpp:2138
int type() const
Returns the type of a matrix element.
_Tp y
y coordinate of the point
Definition types.hpp:202
_Tp x
x coordinate of the point
Definition types.hpp:201
Template class for specifying the size of an image or rectangle.
Definition types.hpp:335
_Tp height
the height
Definition types.hpp:363
_Tp width
the width
Definition types.hpp:362
Template class for short numerical vectors, a partial case of Matx.
Definition matx.hpp:369
Point_< double > Point2d
Definition types.hpp:208
#define CV_16SC1
Definition interface.h:106
unsigned char uchar
Definition interface.h:51
#define CV_8UC3
Definition interface.h:90
@ IMREAD_ANYDEPTH
If set, return 16-bit/32-bit image when the input has the corresponding depth, otherwise convert it t...
Definition imgcodecs.hpp:72
@ IMREAD_LOAD_GDAL
If set, use the gdal driver for loading the image.
Definition imgcodecs.hpp:74
@ IMREAD_COLOR
If set, always convert image to the 3 channel BGR color image.
Definition imgcodecs.hpp:71
CV_EXPORTS_W bool imwrite(const String &filename, InputArray img, const std::vector< int > ¶ms=std::vector< int >())
Saves an image to a specified file.
CV_EXPORTS_W Mat imread(const String &filename, int flags=IMREAD_COLOR)
Loads an image from a file.
int main(int argc, char *argv[])
Definition highgui_qt.cpp:3
GOpaque< Size > size(const GMat &src)
Gets dimensions from Mat.
How to Read Raster Data using GDAL
This demonstration uses the default OpenCV imread function. The primary difference is that in order to force GDAL to load the image, you must use the appropriate flag.
When loading digital elevation models, the actual numeric value of each pixel is essential and cannot be scaled or truncated. For example, with image data a pixel represented as a double with a value of 1 has an equal appearance to a pixel which is represented as an unsigned character with a value of 255. With terrain data, the pixel value represents the elevation in meters. In order to ensure that OpenCV preserves the native value, use the GDAL flag in imread with the ANYDEPTH flag.
If you know beforehand the type of DEM model you are loading, then it may be a safe bet to test the Mat::type() or Mat::depth() using an assert or other mechanism. NASA or DOD specification documents can provide the input types for various elevation models. The major types, SRTM and DTED, are both signed shorts.
Notes
Lat/Lon (Geographic) Coordinates should normally be avoided
The Geographic Coordinate System is a spherical coordinate system, meaning that using them with Cartesian mathematics is technically incorrect. This demo uses them to increase the readability and is accurate enough to make the point. A better coordinate system would be Universal Transverse Mercator.
Finding the corner coordinates
One easy method to find the corner coordinates of an image is to use the command-line tool gdalinfo. For imagery which is ortho-rectified and contains the projection information, you can use the 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
Below is the output of the program. Use the first image as the input. For the DEM model, download the SRTM file located at the USGS here. http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip
Input Image
Heat Map
Heat Map Overlay