1Reading Geospatial Raster files with GDAL {#tutorial_raster_io_gdal} 2========================================= 3 4Geospatial raster data is a heavily used product in Geographic Information Systems and 5Photogrammetry. Raster data typically can represent imagery and Digital Elevation Models (DEM). The 6standard library for loading GIS imagery is the Geographic Data Abstraction Library [(GDAL)](http://www.gdal.org). In this 7example, we will show techniques for loading GIS raster formats using native OpenCV functions. In 8addition, we will show some an example of how OpenCV can use this data for novel and interesting 9purposes. 10 11Goals 12----- 13 14The primary objectives for this tutorial: 15 16- How to use OpenCV [imread](@ref imread) to load satellite imagery. 17- How to use OpenCV [imread](@ref imread) to load SRTM Digital Elevation Models 18- Given the corner coordinates of both the image and DEM, correllate the elevation data to the 19 image to find elevations for each pixel. 20- Show a basic, easy-to-implement example of a terrain heat map. 21- Show a basic use of DEM data coupled with ortho-rectified imagery. 22 23To implement these goals, the following code takes a Digital Elevation Model as well as a GeoTiff 24image of San Francisco as input. The image and DEM data is processed and generates a terrain heat 25map of the image as well as labels areas of the city which would be affected should the water level 26of the bay rise 10, 50, and 100 meters. 27 28Code 29---- 30 31@include cpp/tutorial_code/HighGUI/GDAL_IO/gdal-image.cpp 32 33How to Read Raster Data using GDAL 34---------------------------------- 35 36This demonstration uses the default OpenCV imread function. The primary difference is that in order 37to force GDAL to load the image, you must use the appropriate flag. 38@code{.cpp} 39cv::Mat image = cv::imread( argv[1], cv::IMREAD_LOAD_GDAL ); 40@endcode 41When loading digital elevation models, the actual numeric value of each pixel is essential and 42cannot be scaled or truncated. For example, with image data a pixel represented as a double with a 43value of 1 has an equal appearance to a pixel which is represented as an unsigned character with a 44value of 255. With terrain data, the pixel value represents the elevation in meters. In order to 45ensure that OpenCV preserves the native value, use the GDAL flag in imread with the ANYDEPTH flag. 46@code{.cpp} 47cv::Mat dem = cv::imread( argv[2], cv::IMREAD_LOAD_GDAL | cv::IMREAD_ANYDEPTH ); 48@endcode 49If you know beforehand the type of DEM model you are loading, then it may be a safe bet to test the 50Mat::type() or Mat::depth() using an assert or other mechanism. NASA or DOD specification documents 51can provide the input types for various elevation models. The major types, SRTM and DTED, are both 52signed shorts. 53 54Notes 55----- 56 57### Lat/Lon (Geographic) Coordinates should normally be avoided 58 59The Geographic Coordinate System is a spherical coordinate system, meaning that using them with 60Cartesian mathematics is technically incorrect. This demo uses them to increase the readability and 61is accurate enough to make the point. A better coordinate system would be Universal Transverse 62Mercator. 63 64### Finding the corner coordinates 65 66One easy method to find the corner coordinates of an image is to use the command-line tool gdalinfo. 67For imagery which is ortho-rectified and contains the projection information, you can use the [USGS 68EarthExplorer](http://http://earthexplorer.usgs.gov). 69@code{.bash} 70\f$> gdalinfo N37W123.hgt 71 72 Driver: SRTMHGT/SRTMHGT File Format 73 Files: N37W123.hgt 74 Size is 3601, 3601 75 Coordinate System is: 76 GEOGCS["WGS 84", 77 DATUM["WGS_1984", 78 79 ... more output ... 80 81 Corner Coordinates: 82 Upper Left (-123.0001389, 38.0001389) (123d 0' 0.50"W, 38d 0' 0.50"N) 83 Lower Left (-123.0001389, 36.9998611) (123d 0' 0.50"W, 36d59'59.50"N) 84 Upper Right (-121.9998611, 38.0001389) (121d59'59.50"W, 38d 0' 0.50"N) 85 Lower Right (-121.9998611, 36.9998611) (121d59'59.50"W, 36d59'59.50"N) 86 Center (-122.5000000, 37.5000000) (122d30' 0.00"W, 37d30' 0.00"N) 87 88 ... more output ... 89@endcode 90Results 91------- 92 93Below is the output of the program. Use the first image as the input. For the DEM model, download 94the SRTM file located at the USGS here. 95[<http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip>](http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip) 96 97 98 99 100 101 102
