1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
/*
* gdal_image.cpp -- Load GIS data into OpenCV Containers using the Geospatial Data Abstraction Library
*/
// OpenCV Headers
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/highgui.hpp"
// C++ Standard Libraries
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <vector>
using namespace std;
// define the corner points
// Note that GDAL library can natively determine this
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 dem corners
cv::Point2d dem_bl( -122.0, 38);
cv::Point2d dem_tr( -123.0, 37);
// range of the heat map colors
std::vector<std::pair<cv::Vec3b,double> > color_range;
// List of all function prototypes
cv::Point2d lerp( const cv::Point2d&, const cv::Point2d&, const double& );
cv::Vec3b get_dem_color( const double& );
cv::Point2d world2dem( const cv::Point2d&, const cv::Size&);
cv::Point2d pixel2world( const int&, const int&, const cv::Size& );
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r );
/*
* Linear Interpolation
* p1 - Point 1
* p2 - Point 2
* t - Ratio from Point 1 to Point 2
*/
cv::Point2d lerp( cv::Point2d const& p1, cv::Point2d const& p2, const double& t ){
return cv::Point2d( ((1-t)*p1.x) + (t*p2.x),
((1-t)*p1.y) + (t*p2.y));
}
/*
* Interpolate Colors
*/
template <typename DATATYPE, int N>
cv::Vec<DATATYPE,N> lerp( cv::Vec<DATATYPE,N> const& minColor,
cv::Vec<DATATYPE,N> const& maxColor,
double const& 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;
}
/*
* Compute the dem color
*/
cv::Vec3b get_dem_color( const double& elevation ){
// if the elevation is below the minimum, return the minimum
if( elevation < color_range[0].second ){
return color_range[0].first;
}
// if the elevation is above the maximum, return the maximum
if( elevation > color_range.back().second ){
return color_range.back().first;
}
// otherwise, find the proper starting index
int idx=0;
double t = 0;
for( int x=0; x<(int)(color_range.size()-1); x++ ){
// if the current elevation is below the next item, then 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 - color_range[x].second);
break;
}
}
// interpolate the color
return lerp( color_range[idx].first, color_range[idx+1].first, t);
}
/*
* Given a pixel coordinate and the size of the input image, compute the pixel location
* on the DEM image.
*/
cv::Point2d world2dem( cv::Point2d const& coordinate, const cv::Size& dem_size ){
// relate this to the dem points
// ASSUMING THAT DEM DATA IS ORTHORECTIFIED
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 a pixel coordinate to world coordinates
*/
cv::Point2d pixel2world( const int& x, const int& y, const cv::Size& size ){
// compute the ratio of the pixel location to its dimension
double rx = (double)x / size.width;
double ry = (double)y / size.height;
// compute LERP of each coordinate
cv::Point2d rightSide = lerp(tr, br, ry);
cv::Point2d leftSide = lerp(tl, bl, ry);
// compute the actual Lat/Lon coordinate of the interpolated coordinate
return lerp( leftSide, rightSide, rx );
}
/*
* Add color to a specific pixel color value
*/
void add_color( cv::Vec3b& pix, const uchar& b, const uchar& g, const uchar& r ){
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; }
}
/*
* Main Function
*/
int main( int argc, char* argv[] ){
/*
* Check input arguments
*/
if( argc < 3 ){
cout << "usage: " << argv[0] << " <image_name> <dem_model_name>" << endl;
return -1;
}
// load the image (note that we don't have the projection information. You will
// need to load that yourself or use the full GDAL driver. The values are pre-defined
// at the top of this file
//![load1]
cv::Mat image = cv::imread(argv[1], cv::IMREAD_LOAD_GDAL | cv::IMREAD_COLOR );
//![load1]
//![load2]
// load the dem model
cv::Mat dem = cv::imread(argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH );
//![load2]
// create our output products
cv::Mat output_dem( image.size(), CV_8UC3 );
cv::Mat output_dem_flood( image.size(), CV_8UC3 );
// for sanity sake, make sure GDAL Loads it 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 our output DEM heat map
// Pair format ( Color, elevation ); Push from low to high
// Note: This would be perfect for a configuration file, but is here 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 a minimum elevation
double minElevation = -10;
// iterate over each pixel in the image, computing the dem point
for( int y=0; y<image.rows; y++ ){
for( int x=0; x<image.cols; x++ ){
// convert the pixel coordinate to lat/lon coordinates
cv::Point2d coordinate = pixel2world( x, y, image.size() );
// compute the dem image pixel coordinate from lat/lon
cv::Point2d dem_coordinate = world2dem( coordinate, dem.size() );
// extract the elevation
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;
}
// write the pixel value to the file
output_dem_flood.at<cv::Vec3b>(y,x) = image.at<cv::Vec3b>(y,x);
// compute the color for the heat map output
cv::Vec3b actualColor = get_dem_color(dz);
output_dem.at<cv::Vec3b>(y,x) = actualColor;
// show effect of a 10 meter increase in ocean levels
if( dz < 10 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 90, 0, 0 );
}
// show effect of a 50 meter increase in ocean levels
else if( dz < 50 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 90, 0 );
}
// show effect of a 100 meter increase in ocean levels
else if( dz < 100 ){
add_color( output_dem_flood.at<cv::Vec3b>(y,x), 0, 0, 90 );
}
}}
// print our heat map
cv::imwrite( "heat-map.jpg" , output_dem );
// print the flooding effect image
cv::imwrite( "flooded.jpg", output_dem_flood);
return 0;
}