GEO-registeration of corona images using Landsat 8 satellite images and terrain analysis using DEM
This project looks at the method of geo-referencing corona images of the 1970s with the present day satellite images using the software ERDAS IMAGINE. This was done by plotting control points that are well distributed throughout these images and then resampling it using the bilinear interpolation. Following geo-referencing and resampling, its digital elevation model is developed. It is then analyzed for its terrain characteristics such as slope, aspect and curvature using the software SAGA.
Keywords: Geo-referencing, DEM, slope, aspect, curvature
|DEM||Digital Elevation Model|
|RPC||Resultant Polynomial Coefficient|
|SRTM||Shuttle Radar Topography Mission|
Coordinates of spatial data are referenced to their corresponding positions on the surface of the earth is called geo-referencing of spatial data. This means that a relationship is established between the coordinates on the map or the computer screen and its corresponding latitude and longitude.
Digital Elevation Model is the digital representation of the land surface elevation with respect to any reference datum. They are used to determine terrain attributes such as elevation at any point, slope and aspect.
1.2 Statement of the problems
The regions of Mandi and Kullu being mountainous posed slight difficulties while geo-referencing. Careful observation is required to mark the control points on the same as point of the mountain ridges on the CORONA and satellite images. The slope, aspect and terrain of the DEMs are difficult to read in ERDAS IMAGINE as it is in different shades of grey. So SAGA is used instead.
The objective of this research is to geo-reference the CORONA images of the regions of Mandi and Kullu with satellite images and develop their digital elevation model.
2.1.1 Types of geo-referencing
Geo-referencing can be 1st Order Polynomial, 2nd Order Polynomial, 3rd Order Polynomial and so on. The minimum number of points required for geo-referencing based on the order of polynomial is given by the following formula:
Let ‘t’ be the order of the polynomial, then
Minimum points required =
But we also require a validation point. Thus the total number of points necessary for geo-referencing is given by the following formula:
Total points required=
2.1.2 Types of projection system
A projection is defined as an ordered system of meridians and parallels on a flat surface.
Many types of map projections are being used for map making. They are basically classified into four groups in accordance with the Map Projection Theory or the types of surfaces that are tangent with the globe. The four categories are:
- Equidistant projection
- Equivalent projection
- Azimuthal projection
- Conformal projection
18.104.22.168 Equidistant projection
It is a map projection that maintains scale along one or more lines, or from one or two points to all other points on the map. The distance between the center point of the map and any other point is correct with an equidistant projection.
22.214.171.124 Equivalent projection
The equal area projection retains the relative size of area throughout a map. That is, at any given region in a map, an equal area projection keeps the true size of features.
126.96.36.199 Azimuthal projection
This type of map projection allows a flat sheet to touch with the globe, with the light being cast from certain positions, including the centre of the Earth, opposite to the tangent area, and from infinite distance.
188.8.131.52 Conformal projection
The Stereographic projection has its origin of light on the globe surface opposite to the tangent point. The created curved lines will be defined on more than half of the sphere. The meridians are straight lines adjacent to one another in the central area and become more widely spaced at the map periphery, while the parallels are circles. Shape is maintained in this type of projection, making it applicable for aviation mapping. (Iliffe, J.,2000)
CORONA was created by a small group of CIA, Air Force and private industry experts who were tasked with finding a way to provide broad imagery coverage of the USSR to identify missile launch sites and production facilities. Known to the public as the U.S. Air Force’s Discoverer program, the classified CORONA project operated during the height of the Cold War to collect pictures over the denied areas of eastern Europe and Asia. CORONA also had sister programs: ARGON for mapping imagery and LANYARD, a short-lived program designed for higher-quality imagery. During its operational life, CORONA collected over 800,000 images in response to the national security requirements of the time. On average, individual images covered a geographic area on the Earth's surface of approximately 10x120 miles. On February 24, 1995, Vice President Gore visited CIA Headquarters to announce Executive Order 12951, signed by President Clinton, which released CORONA, ARGON, and LANYARD imagery to the public. (Ruffner, K. C.,1995)
2.1.4 Landsat 8
It was launched on February 11, 2013 and is working till date. Its payload consists of two science instruments: the Operational Load Imager (OLI) and the Thermal Infrared Sensor (TIRS). These two sensors provide seasonal coverage of the global landmass as a spatial resolution of 30m (visible, NIR, SWIR), 100m (thermal) and 15m (panchromatic). Landsat 8 instruments represent an evolutionary advance in technology. OLI is a push-broom sensor with four-mirror telescope and 12-bit quantization. OLI collects data for visible, near infra-red and shortwave infra-red bands as well as panchromatic band. Compared to Landsat 7’s ETM+ bands, OLI provides two new spectral bands, one tailored especially for detecting cirrus clouds and the other for coastal zone observations. TIRS collects data for two more narrow spectral bands in the thermal region formerly covered by one wide spectral band on Landsats 4-7. Landsat 8 is required to return 400 scenes per day (150 more than Landsat 7 is required to capture) and has been providing regularly 550 scenes per day. This increases the probability of capturing cloud-free scenes for the global landmass. Landsat 8 scene-size is 185-km-cross-track-by-180-km-along-track. Cartographic accuracy of 12m or better is required of Landsat 8 data products. (LSDS-1574 Version, 3, 2015)
2.1.5 Dem generation from stereoscopic imagery
Compared to the traditional manual methods that use human operators, automated methods of DEM generation from remote sensing provide efficient, economic and reasonably accurate products covering extended areas of the Earth’s surface. Remote sensing of the Earth’s surface started with photographic film cameras and has been evolving to digital cameras with selective sensing bands, for example, multispectral, thermal, hyperspectral, and radar. As with all stereo 3D reconstruction, two or more images that sense a scene with overlapping areas are required. A stereo image pair can be formed by along-track or across-track arrangement of sensors. Along-track is defined by the forward motion of the satellite along its orbital path, whereas across-track refers to a satellite traveling on different orbits, hence images covering the same area are taken from different orbits. In general, a stereo pair captured in an along-track mission has a shorter time interval between two images than that captured in an across-track mission. Thus variable weather conditions have less effect on along-track stereo pairs than on across-track stereo pairs in these passive imaging scenarios. The distance between two sensors is called the baseline (B), and the nadir distance (vertical distance) from satellite to ground is referred to as the height (H), as illustrated in the figure shown below. The base to height (B/H) ratio is a key parameter in DEM generation from stereoscopic imagery. It is a criterion for choosing an adequate number of stereo pairs from the same or different orbits.
DEM generation from a stereoscopic image pair involves the following processes:
• Pre-processing of image pairs for noise removal: this is a process to mitigate the effects of noise introduced by the image sensors.
• Image matching: this is the process of finding corresponding points in two or more images and is implemented by either area-based or feature-based matching, or a combination of both.
• Triangulation process: image coordinates of matched points from the image pairs are transformed into ground coordinates using the cameras’ interior and exterior parameters. This process involves geometric modeling of the satellite camera system and the ground coordinate system.
• Evaluation of the reconstructed DEM: this process can be achieved by means of ground control points (GCPs), if available. (Pears, N., Liu, Y., & Bunting, P.2012)
2.1.6 RPC file
For the purpose of georeferencing Cartosat-1 Leader file has been used. Ephemeris or the space positional coordinates and the attitude angles, which are a function of time, are extracted from the Leader file. The image coordinates and the corresponding ephemeris coordinates are used to derive the ground coordinates with the help of the physical model (Collinearity equation). The unknown Rational Polynomial Coefficients are obtained by developing a Rational function between the image row, column and the corresponding ground coordinates. From the calculated Rational Polynomial Coefficients, the ground coordinates are derived using the Rational Function. The geodetic height for the corresponding ground coordinates are input from the Shuttle Radar Topographic Mission 90 meter Digital Elevation Model. Then Root Mean Square Error is obtained between the calculated ground coordinates and those obtained by ground survey. (Devika, S., Ramakrishnan, S. S., Mohan, A. M., & Rao, B. S)
The Shuttle Radar Topography Mission, flown on Space Shuttle Endeavour in February 2000 (STS‐99), was a joint project of the National Aeronautics and Space Administration, the National Geospatial‐Intelligence Agency (NGA), and DLR. DLR worked in partnership with ASI. The SRTM objective was to acquire a digital elevation model of all land between about 60° north latitude and 56° south latitude, about 80% of Earth's land surface. In quantitative terms the cartographic products derived from the SRTM data were to be sampled over a grid of 1 arc sec by 1 arc sec (approximately 30 m by 30 m), with linear vertical absolute height error of less than 16 m, linear vertical relative height error of less than 10 m, circular absolute geolocation error of less than 20 m, and circular relative geolocation error of less than 15 m. The relative height error of the X band SRTM data was to be less than 6 m. All quoted errors are at 90% confidence level, consistent with National Map Accuracy Standards (NMAS). These specifications are similar to those of the 30‐m DEMs produced by the U.S. Geological Survey as part of the National Elevation Dataset (NED). NED was produced by photogrammetric reduction of stereo air photographs yielding generally a representation of the elevations of the ground surface even beneath vegetation canopies. The SRTM radars were unable to sense the surface beneath vegetation canopies and so produced elevation measurements from near the top of the canopies.
SRTM employed two synthetic aperture radars, a C band system (5.6 cm, C radar) and an X band system (3.1 cm, X radar). NASA's Jet Propulsion Laboratory (JPL) was responsible for C radar. The DLR with Astrium, its contractor for the X band space segment, was responsible for X radar. The operational goal of C radar was to generate contiguous mapping coverage as called for by the mission objectives. X radar generated data along discrete swaths 50 km wide. These swaths offered nearly contiguous coverage at higher latitudes. X radar was included as an experimental demonstration. As it did not employ ScanSAR, the X band radar had a slightly higher resolution and better signal‐to‐noise ratio (SNR) than the C band system. Thus it could be used as an independent data set to help resolve problems in C radar processing and quality control. (Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., & Seal, D. , 2007)
This is done by plotting required number of control points of the same area on the satellite image and the given corona image, on the ERDAS IMAGINE software, and then resample it. The satellite image is a geo-referenced image while the corona image is not. So by plotting points of the same area on both these images and resampling it will result in the corona image getting resampled.
2.2.2 Digital elevation model
Below is the DEM of the regions of Mandi and Kullu developed in ERDAS IMAGINE. Using this, its slope, aspect and curvature are analyzed.
It represents the rate of change of elevation for each DEM cell. The yellow region represents less slope while the red regions shows high slope.
Aspect identifies the downslope direction of the maximum rate of change in value from each cell to its neighbors. The value of each cell in an aspect dataset indicates the direction the cell's slope faces.
The white region depicts the downslope direction towards the north, the pink region depicts the downslope direction towards the east, the black region depicts the downslope direction towards the south and the blue region depicts the downslope direction towards the west.
184.108.40.206 Planform or plan curvature
Planform curvature is perpendicular to the direction of the maximum slope. A positive value indicates the surface is sidewardly convex at that cell and a negative value indicates the surface is sidewardly concave at that cell. A value of zero indicates the surface is linear. Profile curvature relates to the convergence and divergence of flow across a surface.
220.127.116.11 Profile curvature
Profile curvature is parallel to the direction of the maximum slope. A negative value indicates that the surface is upwardly convex at that cell. A positive profile indicates that the surface is upwardly concave at that cell. A value of zero indicates that the surface is linear. Profile curvature affects the acceleration or deceleration of flow across the surface.
RPC files are required for making DEMs. However, there are no RPC files for CORONA images. Thus, alternate methods of generating DEM should be used.
1. Iliffe, J. (2000). Datums and map projections for remote sensing, GIS, and surveying. CRC Press.
2. Ruffner, K. C. (Ed.). (1995). Corona: America's first satellite program. History Staff, Center for the Study of Intelligence, Central Intelligence Agency.
3. Landsat, U. S. G. S. (2015). 8 (L8) data users handbook. LSDS-1574 Version, 3
4. Pears, N., Liu, Y., & Bunting, P. (2012). 3D imaging, analysis and applications (Vol. 3). New York: Springer.
5. Devika, S., Ramakrishnan, S. S., Mohan, A. M., & Rao, B. S. DEVELOPMENT AND IMPLEMENTATION OF RATIONAL POLYNOMIAL COEFFICIENT ALGORITHMS FOR GEOREFERENCING CARTOSAT-1 DATA.
6. Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., & Seal, D. (2007). The shuttle radar topography mission. Reviews of geophysics, 45(2).
I would like to thank my professor Dr. Dericks P. Shukla for giving me the opportunity to learn and understand about remote sensing and for guiding me. I would also like to express my sincere gratitude to Mr. Sharad and Mr. Niraj for helping me throughout my project.
I am extremely grateful to the Indian Academy of Sciences and the Indian Institute of Technology Mandi for providing me with the summer research fellowship.
Most importantly, I feel extremely blessed to have had the support of my parents and sister till the very end.