• CorrGeom: Geometric correction of raster and vector files

Access this help text as a web page: CorrGeom

Presentation and options

Introduction

Conventional paper maps give us an enormous amount of information which can be incorporated in to a Geographic Information System (GIS). Although more and more maps are now generated digitally there are still many cases in which digital files are not available, either because we are dealing with documents of a certain age, or because they have been created using conventional analogue techniques, or simply because we don't have access to the original digital information.

In these cases, although digitisation using a digitising table may be a possibility it requires a laborious effort and a loss of detail compared to the possibility of scanning the digital document and proceeding to its digitisation on the computer screen.

MiraMon allows documents from practically any scanner (in formats BMP, TIFF, JPG, etc) to be transformed to IMG or JPG formats (if they were not already in one of these formats). Once in these formats it is possible to digitise over them the contour lines, roads, land use polygons, etc and convert them in to vector layers. As the scanned documents are still not georeferenced when we digitise directly over them the resultant vector coordinates will be in some arbitrary units (pixels) of the scanning system and not in the map units of the desired georeferenced coordinate system (for example, in meters in the UTM system). Although it is possible to proceed in this way and to georeference the resultant vector document after digitising using the MiraMon CorrGeom program, it is advisable to georeference the raster first and then to proceed to its digitisation. This will also allow map sheets to be scanned separately and then precisely mosaiced together as well as permitting them, once georeferenced, to be displayed with any other available digital vector information layer overlain for reference. Georeferencing of scanned raster files can also be performed using the MiraMon CorrGeom program. So, in summary, there are two alternatives for introducing vector layers into GIS from scanned paper documents using the CorrGeom MiraMon program:

• Scan -> Digitise -> Georeference the vectors
• Scan -> Georeference the raster ->Digitise (RECOMMENDED)

In both cases the recommended method is the first order polynomial fit. In the very unlikely case, given the quality of modern scanners, that the scanned images show quadratic distortions, use the second order polynomial fit.

There exist a number of other interesting sources of information for incorporation in our GIS: such as raster layers like remotely sensed images and aerial photography. These data can be obtained directly in digital format (this is the case for the majority of remotely sensed images) or, alternatively, they may be available in paper or photographic form (the usual case for conventional aerial photography) which can be scanned to convert them into digital images.

Remotely sensed images present a series of geometric distortions produced by the rotation and curvature of the Earth, the irregularities of the platform's orbit, etc. The MiraMon CorrGeom program allows such images to be corrected to remove these distortions using conventional techniques based on first and second order polynomials. In the case of aerial photographs and remotely sensed images with a high level of detail (50m or smaller pixel size) it is advisable to perform corrections taking into account the effect of the relief of the land surface. CorrGeom offers the methods "polynomial fit with heights" for satellite images with small pixels (Landsat TM, SPOT XS and P, IRS LISS, etc) and "colinearity equations fit" for aerial photographs which applies fundamental photogrammetric equations. The images that have been corrected using information about the relief contained in a Digital Elevation Model (DEM) are known as orthoimages (or orthophotos in the case of aerial photos).

The CorrGeom program therefore allows the correction of these distortions (due to the platform's orbit, the relief, etc) and adapts the images to a known cartographic projection system, such as UTM. Once in the appropriate projection system the images can be overlain with other images and vector layers referenced to the same projection system.

Whether correcting raster or vector data it is usual to have a set of ground control points (GCP) which indicate a number of coordinate points in both the uncorrected original reference system and in the final, corrected, reference system. The GCPs can be input with a text editor or digitised as points in MiraMon and converted into GCPs with VECCOR (see below).
In the case of images acquired with a system for which the position and attitude of the sensor are known with precision (GPS+INS) it is possible to do the correction without GCPs.

To increase the statistical reliability of the geometric correction process CorrGeom allows a control points file to contain two subsets of points, one that is used to determine the transformation equations for the fit and the other, which we call the test points, that is used to estimate the error of the fit based on an independent set of points. For both sets of points CorrGeom provides the RMS error in X, Y and globally for each point as well as for the set of control points. Logically, if only one set of points is given then the RMS errors are only given for this set.

If the user wishes to change the projection system, whether of raster or vector datasets, once they have been georeferenced with MiraMon (or if they were imported into MiraMon with georeferences) then the MiraMon program CanviPrj can be used. Do not confuse the functions of CorrGeom that allows non-georeferenced layers to be georeferenced and CanviPrj that changes the projection of layers for which the georeference is already known.

NOTES:

1. When we speak of a projection system we are referring to UTM, Lambert, etc, whilst when we speak of a reference system we are referring to, for example, UTM-31N-UB/ICC (ie. the UTM with huso 31, Northern hemisphere, European ED50 datum with the UB/ICC adjusted parameters), or Goode_Homolosine-WGS84 (Goode Homolosine WGS84 with the WGS 1984 datum), etc. (see the help on Geodesia for more information)
2. From now on we refer to images and vectors that have been adapted to some reference system using CorrGeom as "corrected".

CorrGeom

The CorrGeom application allows the geometric correction of rasters (IMG and JPG: satellite images aerial photographs, scanned maps...) or vector layers (VEC, PNT, ARC and their derivatives POL and NOD) using a series of n known control points:

• Coordinates in the original reference system ("uncorrected system"), eg. in pixel units, typically with the origin at (0,0) in the bottom left corner in rasters, or eg. in cm in a file coming from the digitisation of a vector paper map created on a digitising tablet and calibrated in centimetres on the map sheet. We will call these units original reference system coordinates: X_SistOri, Y_SistOri.
• Coordinates of the chosen reference system ("corrected" or "map" system, for example UTM-31N). These coordinates are determined manually with the help of map sheets, on the screen with other digital maps, with GPS in the field, etc. We will call these coordinates reference system coordinates : X_SistDest, Y_SistDest [, Z_SistDest].

In addition, the program allows to apply affine transformations (translations, rotations, tilts, scaling and mirrors) and perspective to both rasters and vectors.

All the bands of a multiband raster (typically a 24 bit RGB or a multispectral satellite image) can be geometrically corrected simultaneosuly by using the /MULTIBANDA parameter or by activating the corresponding button. All the bands will then be corrected in a single operation whilst maintaining a multiband raster as the output. In this case we recommend editing the control points using the band with the highest resolution. If the correction requires the use of a digital elevation model (DEM) and the different bands have different pixel sizes it is possible to use a multiband DEM with as many bands as there are different pixel sizes by activating the parameter /MDE_MULTIBANDA.

The correction can be performed using various methods:

First degree polynomial fit: This method does not consider height information (elevation). Typically it is used to correct vector and raster data with little distortion, such as scanned images and maps that were previously corrected (like orthophotos or topographic maps) or satellite images in areas with little marked relief, limited extent, and a relatively large pixel size in relation to the altitude of the satellite. If the objective is to obtain cartography with high geometric quality then the application of these methods to scanned aerial photographs or relatively detailed and extensive satellite images (like SPOT or Landsat-TM) is only appropriate in regions with flat relief.

Polynomial fit with Z for the columns and first order polynomial without Z for the rows: With this method the fit of the columns is made using a polynomial that simultaneously takes into account X, Y and Z. As for the rows, the fit can be made in three different ways: considering that there is no pitching movement of the sensor during the acquisition of the image, considering that the pitch is constant or considering that the pitch varies linearly. For more details see:
Palà, V., Pons, X. (1995) Incorporation of relief into geometric corrections based on polynomials. Photogrammetric Engineering & Remote Sensing, 61(7):935-944.

As for 3 but with a constant pitch and Z for the rows: With this method the fit of the columns is made using a polynomial that simultaneously takes into account X, Y and Z. As for the rows, the fit can be made in three different ways: considering that there is no pitching movement of the sensor during the acquisition of the image, considering that the pitch is constant or considering that the pitch varies linearly. For more details see:

Palà, V., Pons, X. (1995) Incorporation of relief into geometric corrections based on polynomials. Photogrammetric Engineering & Remote Sensing, 61(7):935-944.

As for 3 but with a variable pitch and Z for the rows: With this method the fit of the columns is made using a polynomial that simultaneously takes into account X, Y and Z. As for the rows, the fit can be made in three different ways: considering that there is no pitching movement of the sensor during the acquisition of the image, considering that the pitch is constant or considering that the pitch varies linearly. For more details see:

Palà, V., Pons, X. (1995) Incorporation of relief into geometric corrections based on polynomials. Photogrammetric Engineering & Remote Sensing, 61(7):935-944.

Colinearity equations fit: The parameters used in the colinearity equations are: position (Xc, Yc, Zc) and attitude (omega, fi, kappa) of the camera or sensor, the camera focal length (f) and the coordinates of the principal point (x_pp, y_pp) of the image. The user may decide which of these parameters to estimate by fitting and which are fixed. For those which are to be estimated it is useful to give some approximate values which will help the program converge more rapidly. In some cases, if the initial parameter is completely unknown or the approximation is very poor, the algorithm may not converge or unsatisfactory solutions may be produced. If all the parameters are known and fixed then the transformation is performed without using control points and using the test points just to evaluate the error of the fit. These are the initial values of the flight parameters in the first iteration:
• Roll: omega = 0
• Pitch: fi = 0
• Yaw: kappa = Various values are tried and the one that gives the smallest RMS in the first iterations is used.
• Xc = Average of the X coordinates of the control points (in the reference system).
• Yc = Average of the Y coordinates of the control points (in the reference system).
• ZZc = Average Zc + flight altitude, approx. 2000-3000 m
When the principal point (x_pp, y_pp) has to be fit the program takes as initial values:
• x_pp = Average of the x coordinates (in original units) of the control points or the middle of the raster.
• y_pp = Average of the y coordinates (in original units) of the control points or the middle of the raster
THowever, sometimes these averages can be quite distant from the real principal point and in these cases the program may not converge (eg if only a fragment from the edge of an aerial photograph is scanned). The program allows the user to change the value for another estimated by eye at the intersection of the lines joining the fiducial points. Remember to use always the same measurement system, consistent with the original reference units.
For fitting the focal length it is necessary to choose some initial value close to the real value, typically between 15 and 35 mm for digital cameras and between 50 and 250 mm for analogical cameras.
In addition it is necessary to know the number of points per inch (ppp, ppi or dpi) of the scanned image (1 inch = 25.4 mm). For analogue photos common values are between 400 and 2540 dpi (exceptionally 3175). For digital photos refer to the manufacturer's technical specifications (for example the AA497-AMDC 28.0 camera with 2024 columns x 2041 rows and a sensor size of 18x18 mm gives a resolution of 2811-2835 dpi). To convert from scanning resolutions given in dots per centimeter to dpi, multiply the value by 2.54. To convert from scanning resolutions given in micrometers (µm) to dpi, divide 25400 by the value.

Geometric Transformations: These types of transformations follow a different philosophy from the previous ones: instead of adjusting the transformation using control points, we want to specify the transformation parameters ourselves because we know them a priori. For instance, we may have scanned a document rotated 90º to fit better into the scanner once rotated; in this case, it will be easier to indicate that we want to apply a 90º rotation to return it to its original position and view it in a "natural" way (perhaps before placing control points for fine-tuning the geometry). A second example would be having an old vector dataset from software that stored coordinates with single precision (such as PC-Arc/Info) and, for this reason, truncated the most significant digit when it was constant for the entire layer (for example, in Catalonia in UTM-31N, it was common to omit the 4 from the Y coordinate, so a coordinate like 4619254.734 was written as 619254.734, maintaining decimetric precision); in this case, we can recover the original coordinates by specifying a translation transformation of magnitude 4000000. Finally, a third example would be having a dataset where the horizontal reference system units are in km, and we want to convert them to meters, which requires applying a scaling factor of 1000.
The types of operations offered are translations, scalings, shears, rotations, perspectives, and mirrors. For affine transformations (all except perspective), a transformation origin can be specified, and for perspective, an envelope over which to act can be specified. The possibility of choosing this is provided so that if we want to modify a vector layer over an already modified raster, we can use the raster's envelope to perfectly overlay the transformed vector layer onto the modified raster.
It is possible to perform a sequence of transformations and/or a perspective one after the other on a raster or vector, with the sole restriction that only one perspective can be used.

Adjustment based on reference points: In this case, the geometric correction is made from two rasters that, when paired, provide the origin position of each pixel of the destination image (therefore, they typically have the same dimensions as this one) and together they form, therefore, a mesh of reference.
The procedure is based on the following: For each pixel center of the desired output raster (rectified, corrected, georeferenced, depending on what the coordinates are used for), the program searches for the values of this position in the reference rasters (one for the axis of the X and the other by the Y axis), which ends up providing the specific pixel in the image to be corrected that is closest to the coordinates of the reference rasters. The value in the source raster will be assigned to this output raster pixel when the rescan method is nearest neighbor, or it will be interpolated from neighbors in the other available rescan methods. By default, the output file will have the same envelope, pixel side, and reference system as the XY reference rasters.

COR files contain the ground control points and, optionally, test points, in the following format:

 Number_of_points [Comments]
 X_OriSyst1 Y_OriSyst1 X_DestSyst1 Y_DestSyst1 [Z_DestSyst1] [Comments1]
     :          :          :           :            :            :

• Z coordinate Z_DestSyst1 are not required in the first two options of the program.
• Test points (optional) will be written, with the same format, below this sequence. The user must write at the start of the relevant block the number of test points.
• The file may end with comments of indefinite length, which the program does not preserve and which are replaced by the results of the adjustment.
• The "Original System" coordinates must be Cartesian, i.e., the coordinates must increase positively rightwards in x and upwards in y. The coordinate 0,0 can be either at any point of the original file or even outside it.

This format can be generated with a text editor or from a file of points digitised over the original file. This file can be easily transformed to COR format using the VECCOR program.

For raster images and for program options that do not require input of either the resolution or the bounding rectangle of the output file, these parameters are taken from the DEM documentation file. The maximum envelope values of the output file are adapted as a function of the value of the output resolution (res) using the following expression:

    Xmax = Xmin + res * rounding_by_ceiling (Xmax-Xmin/res)
    Ymax = Ymin + res * rounding_by_ceiling (Ymax-Ymin/res)

The number of rows and columns of the output image can be determined from:

    Ncolumns = rounding_by_ceiling (Xmin-Xmin/res)
    Nrows = rounding_by_ceiling (Ymax-Ymax/res)

The background or NoData value will be written in the corrected filewhen the uncorrected file has no data for the output cells in that position.The user may choose the NODATA value or allow the program to set this. Bydefault it will be the same as the input, if this exists or, if there are noNODATA values in the original data, the default is 255 for byte/byte-RELimages or the smallest allowable negative value for all other formats. If theimage to be corrected is byte/byte-REL and has no defined NODATA value, and ifwhen it is corrected NODATA values are generated then by default the originalvalid values that in the output correspond to NODATA will be saturated. In 24bit images reduced to 8 bits with an optimised palette this can lead to anincorrect visualisation so consequently an option is given for the program toreclassify these values to the closest colour index so that the displayremains practically unchanged.

The program supports digital elevation models of type BYTE, INTEGER (short) and REAL, compressed and uncompressed.

Dialog box of the application

CorrGeom dialog box

Graphic examples

Detail of the placement of a control point at the UTM-06N coordinate (453000,7194000), North America datum 1983, on a 1:25000 USGS topographic map in Alaska Fairbanks. The careful placement of a point at each of the 4 ends of the map allows to obtain, georeferenced, the cartographic document with an RMS that has an adjustment of 28 cm with the scanned document (in reality, since the map is at a scale of 1:25000 the real RMS on the ground must be about 5 m).

Adjustment points located on a SPOT image ((c) SPOT-Image)	General aspect of the correction of the previous SPOT image on an ICGC orthophoto (areas with NoData are shown transparent).
Detail of the correction of the previous SPOT image compared to an image from the ICGC. The reference grid has a side of 1 km.

Adjustment (A) and test (T) points located in a frame from the Cartographic and Geological Institute of Catalonia (ICGC).	Detail of the correction of the previous frame, converted into an orthophoto (left) compared with an orthophoto from the ICGC (right) in order to appreciate the very good geometric fit achieved. Each pixel in the image corresponds to 50 cm.

Syntax examples

For raster images:

CorrGeom 1 MAP MA2 LANDSAT.img  LND 0 450000 480500 4724500 4750000 30 /NODATA=0
CorrGeom 2 MAP NUL LANDSAT.img  LND 1  X  X  X   X  30
CorrGeom 4 MAP MA2 LANDSAT.img  ORTOLND 1 MDE30 
CorrGeom 6 MAP MA2 FOTOSCAN.img ORTOFOT 1 MDE05 400 20 153 1475 1512

For vector files

CorrGeom 1 MUNIPC_1 MUNIPC_2 MUNI_1.pnt MUNI_COR
CorrGeom 2 MUNIPC_1 NUL  MUNI_1.arc MUNI_COR
CorrGeom 1 MUNIPC_1 MUNIPC_2 MUNI_1.vec MUNI_COR 0

For geometric transformations (afines and perspective):

6.5 units translation in positive vertical direction:

CorrGeom 7 LANDSAT.img  LND 0 T 0 6.5

7 units translation in negative horizontal direction:

CorrGeom 7 LANDSAT.img  LND 0 T -7 0

45 degrees clockwise rotation:

CorrGeom 7 LANDSAT.img  LND 0 R 45

45 degrees anticlockwise rotation:

CorrGeom 7 LANDSAT.img  LND 0 R -45

Scaling ratio 2 (in both directions: vertical and horizontal) plus 90º clockwise rotation:

CorrGeom 7 LANDSAT.img  LND 0 E 2 2 R 90

Slight inclination in horizontal right direction:

CorrGeom 7 LANDSAT.img  LND 0 I x 0.1

50% inclination in horizontal left direction:

CorrGeom 7 LANDSAT.img  LND 0 I x -1

Perspective of points in a raster which has a perspective (of 0.5).
Supose that the envelope of the raster is (XMIN, YMIN)=(0,10), (XMAX, YMAX)=(4500, 6000)

CorrGeom 7 MUNI_1.pnt MUNI_COR P 0.5 /XMIN=0 /XMAX=4500 /YMIN=10 /YMAX=6000 Notice that if we don't apply the optative parameters, these points won't adapt to the modified raster.

Rotation of points in a raster which has been also rotated (45º).

CorrGeom 7 MUNI_1.pnt MUNI_COR R 45 /XORI=0 /YORI=10

Notice that if we don't apply the optative parameters, these points won't adapt to the modified raster.

Scaling from cm to m and after translate 4000000 m in Y's. The arguments of affine transformations are transformed for the program itself, so, we can think this operation as: Scaling from cm to m. Traslate 400000000 cms in Y's.

CorrGeom 7 Fitxer.arc Fitxer_sor.arc E 100 100 T 0 400000000

For adjustment based on reference points:

CorrGeom 8 Fitx_in.img Fitx_out.img Fitx_X.img Fitx_Y.img 0 x x x x
CorrGeom 8 Fitx_multi_in.img Directori Fitx_X.img Fitx_Y.img 0 x x x x /MULTIBANDA /PREFIX=apr_

CorrGeom  1/2  cor cor_s  f_i.img  f_s int xm  xM  ym  yM   res
CorrGeom 3/4/5 cor cor_s  f_i.img  f_s int MDE
CorrGeom   6   cor cor_s  f_i.img  f_s int MDE dpi iter op_posicio op_actitud op_focal op_ppral
CorrGeom   7   f_i.img  f_s int [T Tx Ty][E Ex Ey][I Id Ir][R g][P Pr][M Md][...]

CorrGeom 1/2 cor cor_s f_i.pnt f_s
CorrGeom 1/2 cor cor_s f_i.arc f_s
CorrGeom   7 f_i.pnt  f_s [T Tx Ty][E Ex Ey][I Id Ir][R g][P Pr][M Md][...]
CorrGeom   7 f_i.arc  f_s [T Tx Ty][E Ex Ey][I Id Ir][R g][P Pr][M Md][...]

CorrGeom 1/2 cor cor_s f_i.vec f_s [n_decimals]
CorrGeom   7 f_i.vec  f_s [T Tx Ty][E Ex Ey][I Id Ir][R g][P Pr][M Md][...]