Back to Notes_on_Map_Geometry_and_Projections

__What are Coordinates?__

Coordinates or better Co-Ordinates is a term used to express
that two or more information stream relate to each other.

*Strictly speaking, we need to include a single information
stream, a single Ordinate for completeness. Take a calendar
for example. A calendar is a linear stream of days sequentially
arranged one after the other. Graphically we can mark each day with
it's date (day/month/year) or the name of the month or week or day
along a straight line, the ordinate line. We can also mark special
days for the event falling on the day along this line (i.e. your
birthday).*

True Coordinates of 'two dimensional (2D)', 'three dimensional (3D)' or higher 'n-dimensional' character relate as many information stream to each other as the dimensional number expresses.

A simple example of a 2D coordinate system is a language dictionary in form of a 'Look-up Table': alphabetically ordered words of one lanuage are related to the equivalent words of another language. (Note: this kind of a coordinate system can also be interpreted as a simple (co)-ordinate transformation: jumping from one word in one language to the respective word in the other does not change the meaning, it just changes the way how this meaning is expressed.

Examples of 2D and 3D coordinates systems more alingned with GIS
applications are those used to express positions in form of
*logitude* and *latitude, plane map coordinates X and Y, 3D
terrain coordinates X, Y and Z, screen coordinates x and y*,
etc.

The position coordinate examples named above are of two types:

- Polar coordinates: involving one or two angles and a radial distance
- Rectangular cartesian coordinates

In GIS applications we sometimes also encounter a third type: non-rectangular cartesian coordinates (more about this in the section about Coordinate Transformations).

**Polar Coordinates** [1]

(Source: en.wikipedia.org/wiki/Polar_coordinate_system#Position_and_navigation)

In __ 2D polar coordinates systems__ the polar angle
and the radial distance are defined on a surface (plane or curved).
The radial angle may be defined clockwise (normal in GIS) or
anti-clockwise (as shown above).

Polar coordinates for points in the landscape are most often the result of conventional surveying methods (in a horizontal plane): The polar axis L is the direction to North (or can be related to it), the radial distances are horizontal distances (or can be related to those). Angles are measured in the field using Theodolites (or computed from map coordinates), radial distances are measure in the field using steep tapes or Electronic Distance Measurement instruments (EDM) (or can be computed from map coordinates).

Conventional polar surveying results can be transformed into rectangular cartesian coordinates using trigonometric functions: Assume L points to geographic North and angles are measured clockwise, then

X(point) = X(O) + r x sin(angle) ==> pointing eastward

Y(point) = Y(O) + r x cos(angle) ==> pointing northward

__ 3D Polar Coordinate systems__ are in GIS applications
most often encountered in form of 'Geographic Coordinates' of
latitude and longitude (the 3rd dimension is assumed to be the
radial distance between a points location on the ellipsoid and its
gravitational vertical's intersection with the polar axis). Using a
spherical earth model this intersection point would be Earth's
center, an ellipsoidal model complicates this.

Logitudes are defined as angle measured in Earth's equator plane starting at the 'Prime Meridian' (the meridian passing through Greenwich observatory) and the meridian passing through the referenced point as 0 to 180 degree East or West. Often eastern longitude are shown as positive angles and western longitudes as negative values.

Latitudes are defined as angle measured in the meridian plane of the meridian passing through the referenced point between the equator plane and the gravitational 'normal' (vertical) of this point from 0 degree for a point at the equator to 90 degree at the North Pole and -90 degree at the South Pole. Latitudes are often used as un-signed angles but with an appended letter 'N' or 'S' indication north or south of the equator.

Example: a place near Brisbane, Australia, may be referenced as

153E 30S or

153 -30

units are in degree, minutes and seconds or in decimal degree, the '+' normally is not shown.

**Rectangular Coordinates**

For more information see [2] or similar.

For mapping and GIS applications the following is preferred:

A. __2D coordinates (x,y)__

- Plane x, y coordinates (dominant case):
- x and y axis are straight ordinate axis spanning open a plane surface

- (curved) surface x, y coordinates:
- x and y axis are straightest ordinate axis on a curved surface (ie. on a globe)

- x and y axis intersect in a point called the 'Origin' of the coordinate system
- x and y axis intersect at right angle
- distances along each axis are referenced against this origin : x=0 and y=0
- distance scale is set either as linear (metric; preferred case) or as a non-linear function of the distance itself (e.g. logarithmic)
- preferred axis arrangment: When drawing x and y axis on a sheet of paper, the x axis points with positive x values to the right and the y axis with positive y values upward (away from you)
- (see above sketch on the left)

Examples:

- UTM map coordinates:
- Axis Names: x-axis =
**Easting**; y-axis =**Northing**when using False Origin - True Origin: intersection between Central meridian and Equator
- False Origin: produced by adding a constant value to x and y referenced against True Origin; often Easting = x + 500000 m and Northing = y + 1000000 m

- Axis Names: x-axis =
- Scaled UTM map sheet coordinates: Paper and digital maps are
derived from UTM coordinates of a UTM mapping Zone. These zones
extend from high northern to high southern latitude, nearly 20000
km on Earth. Systematically arranged square or rectangular sections
are show on paper maps at a desired scale (e.g. 1 : 100000, 1
cm on the map represents 100000 cm on Earth. The latter is the same
as 1 km on Earth; a 50 cm by 50 cm map shows an area of 50 km by 50
km on Earth). -- If you work in imperial units, replace 'cm' with
'"' --.

Taking measurements on such a map one would use an appropriate graded scale ruler or just a mm scale. Scaling coordinates of a map will refer to the bottom left corner of the map as local map origin (LMO), reading point coordinates as x' and y' for examples in units of mm.

UTM Zone coordinates then can be found as- UTM-x = UTM-x(LMO) + x' * RN with RN = Representative scale
Number (e.g. 100000/unit convertion factor; if x' is in mm and
UTM-x is in m, this factor is 1000)

UTM-x = UTM-x(LMO) + x' * 1000 and UTM-y = UTM-y(LMO) + y' * 1000

- UTM-x = UTM-x(LMO) + x' * RN with RN = Representative scale
Number (e.g. 100000/unit convertion factor; if x' is in mm and
UTM-x is in m, this factor is 1000)
- Computer Screen coordinates:
- Origin at top most, left most screen pixel (x=0; y=0)
- x-pixel count from origin to the right pixel by pixel
- y-pixel count from origin downward (
**!**) scan line by scan line

Opposite to normal y axis direction in mapping and GIS

- Program Window coordinates:
- Origin at programmetrical fixed screen position, often x > 0 and y > 0, for example at screen pixel 100 on screen scan line 150
- Otherwise the same as for Computer Screen coordinates.

- Special cases of plane coordinate use:
- Rater and/or Grid data:
- Raster Images:

Raster Images come in various flavours and formats, like BMP, JPG, PNG, TIF, etc.

Files of Raster Images generally contain two information types (a) a raster and file format description and (b) the image data itself.

The image data in general is organized very similar to that of a computer screen image; how to interpret this data set is specified in the raster and format description. In most raster images, the pixel organisation is assumed as pixel by pixel forming image scan lines and line by line to form the image. No explicite x-pixel and y--pixel positions are given. It is assumed that the sequence of the image data in the file follows this positional sequence and can be re-established on the fly.

For detailed descriptions of Raster Image data formats, see the internet for the specific file type. - Digital Evelation Models (DEM):

DEMs are again very similar in format to Raster Images. DEM files usually come in one of two forms (a) ASCII or (b) binary. They often carry two sections analog to Raster Images, (1) a short descriptive section and (2) the data itself.

The short descriptive section carries information like ground resolution cell (pixel) size, number of cells per scan line, number of scan lines, top left point's map coordinates (e.g. UTM) and a default missing data value in cells (often -9999). No other positional information is given otherwise.

- Raster Images:

- Rater and/or Grid data:
- Photo coordinates used in photogrammetric restitution of
(aerial, terrestrial) photos [3]
Photogrammetric methods are used to reconstruct terrain shape from
stereo photographs. These methods require a strict geonetric
definition of several related coordinate systems. Starting with the
geometry of the camera and the exposed image, coordinate
transformations lead to 3D Earth or Map bound results.

Below a brief outline sketch (terrestrial case) for a single photo situtation:

and a sketch showing the principle of stereo photogrammetry (intersection of Vectors in space).

(Source for both: --Gngdowid 02:33, 28 December 2009 (UTC)) For aerial photos imagine the camera(s) pointing downward (paralell to earth-bound 'Z' axis). The photo coordinate system (one each per photo) have their origin in the socalled 'Principle Point' of the camera system, the x' axis points in direction of aircraft flight and the y' axis completes a right hand plane coordinate system.

B. __3D coordinates (x,y,z)__

Current GIS software is not designed to handle 3D positional coordinates. However, software developers are making attempts to include the 3rd dimension (Height) through using them as special attribute and built-in or add-on (plugin) functions.

In naturally 3D applicationa as in geology or mining, where strata and underground mining levels need to be combined into a single map data set, multiple 3D layers are a crutch to an end.

More advanced GIS software packages like ArcGIS by ESRI use 3D DEM
information to produce 3D visualisation, but are still by nature 2D
GIS applications.

In this sense, a DEM is a 2D information set with the 3rd dimension attached as an attribute.

__Coordinate Transformations__

GIS software like MapWindow, ArcGIS, FGis, GRASS, etc. are coded around several action areas:

- Operating Environment: GUI (Graphical User Interface), Data I/O (basic file access and screen display) Program language support (libraries)
- Database Management System (DBMS): SQL, MYSQL, Postgre, etc. for programmed and user defined data access and interogation
- GIS specific
**Transformations of positional data** - Generation of application specific and/or map specific output

In this section **Transformations of positional data** is of
particular interest.

The 'Simplified Data Flow' of handling positional information is
built into GIS software and can remain hidden from the user. The
only exception is where to find and interpret information relating
to the first two steps.

You find this information in in files of the type *.prj, *.mwprj,
etc.

Below examples of the relevant information of a MapWindows Project file (*.mwprj) and a shapefile (*.prj):

From a MapWindow Project file (*.mwprj):

<Mapwin name="MapWindow+GIS++-+Gatton_2008*"
type="projectfile" version="4.7.3"
ConfigurationPath="..\..\..\..\..\..\Documents and
Settings\gdowideit\Application Data\MapWindow\mapwindow.mwcfg"
**ProjectProjection="+proj=utm +zone=56 +south +ellps=aust_SA
+units=m +no_defs " MapUnits="Meters"**
StatusBarAlternateCoordsNumDecimals="3"
StatusBarCoordsNumDecimals="3"
StatusBarAlternateCoordsUseCommas="True"
StatusBarCoordsUseCommas="True" ShowFloatingScaleBar="False"
FloatingScaleBarPosition="LowerRight" FloatingScaleBarUnit=""
FloatingScaleBarForecolor="-16777216"
FloatingScaleBarBackcolor="-1" MapResizeBehavior="0"
ShowStatusBarCoords_Projected="True"
ShowStatusBarCoords_Alternate="Kilometers"
SaveShap>Settings="False" ViewBackColor_UseDefault="True"
ViewBackColor="16777215">

The projection type from geographic to UTM coordinates is shown
above in bolded letters.

And here a section of a layers *.prj file:

PROJCS["**UTM Zone 56, Southern
Hemisphere**",GEOGCS["**Australian Natl** & S. Amer.
1969",DATUM["D_unknown",**SPHEROID["aust_SA",6378160.0,298.25]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",500000.0],PARAMETER["false_northing",10000000.0],PARAMETER["central_meridian",153.0],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0.0],UNIT["Meter",1.0]]**

Indicating as the above *.mwprj file that this shapefile is based
on a UTM projection referenced against the Australian National
Spheroid from 1966 (now out dated).

As long as the information input is from shapefiles with project
(*.prj) information known to the GIS software package, the
transformation from inupt coordinate system to that specified as
the project coordinate system is performed 'on the fly' in most GIS
packages and hidden to the user.

However, in cases were for example images from aerial photography
or satellite images are to be added, the projection information is
unknown and will have to be determined by user action, called image
registration or 'Georeferencing of images'. The result is a set of
'Transformation Parameters' but no change to the image itself. The
change, display at correct scale, position and orientation is
achieved in GIS software by displaying the image pixels according
to what is specified in the Transformation Parameter set.

Below such a set of Transformation Parameters as can be found in a secondary file alongside the original image file:

1.10401683181001

0

0

-1.10401683180987

432638.106766178

6953073.61784677

Lets refer to these values by names:

**c** = 1.10401683181001

**d** = 0

**e** = 0

**f** = -1.10401683180987

**a** = 432638.106766178

**b** = 6953073.61784677

The transformation then is defined as

**Easting = x(UTM) = a + c x' + d y'** and

**Northing = y(UTM) = b + e x' + f y'**

where x' and y' are image pixel count positions.

The result is that after this transformation pixel positions are
know in map coordinates (here UTM - southern hemisphere). The
transformation is usually done on the fly while displaying the
image (process is hidden from the GIS user). At the same time
(subsequently for each pixel) transformations into program display
window coordinates takes place to finally show the pixel(s), the
image on screen.

**Interpreting the two transformation equations:**

**X = a + c x' + d y'** and

**Y = b + e x' + f y'**

where x' and y' are image pixel count positions and X and Y are the
resulting map coordinates (of the selected map projection
type).

These two equations apply actions of coordinate shifts (a and b),
rotation (alpha_x and alpha_y) and scaling (sx ans sy) (c, d, e and
f).

The total number of parameters used here is 6 (six), two shifts, two rotations and two scales, hence the Name 'Six-Parameter Transformation'.

Another name for the unmodified Six-Parameter Transformation is
'Affine Transformation', as it caters for different scales along
the two coordinate axis and different amounts of axis rotation.

MapWindow with its Georeferencing plugin currently uses a modified approach where

**c = f**- and
**e = -d**

This approach reduces the parameters from six to four ==> 'Four Parameter Transformation'. This tarsformation is also known under the name 'Similarity Transformation'.

ArcGIS with its Georeferencing function allows for both, Four and
Six Paramter Transformations plus for non-linear higher order
polynomial terms to simulate geo-referencing by 'Rubber-Sheeting'
(catering for non-linear image or map distortions).

With this let's interpret the parameter example given above:

**c** = 1.10401683181001

**d** = 0

**e** = 0

**f** = -1.10401683180987

**a** = 432638.106766178

**b** = 6953073.61784677

The first thing to note is that **c = f** is violated!

Why? Each of these contribute to different map coordinates using also different image coordinates as multiplier. The transformation equations given above assume that the arrangement of the coordinate axis to each other is the same. In our case of registering a digital image (y' increments downward) onto a map (Y increments upward) the violation is exaplained.

The second equation should read **Northing = y(UTM) = b + e x'
+ (-f) y'** because of the direction reversal in
y-coordinates.

The condition of **d = e** is not violated. Minus Zero is
still Zero ==> nothing.

In order to learn about scales and rotation amounts we need to
reveal what is hidden inside the four parameters c to f:

- c = sx * cos(alpha_x)
- d = sy * sin(alpha_y)
- e = sx * sin(alpha_x)
- f = sy * cos(alpha_y)

Inspect the above with trigonometric identities in mind, for
example

tan(angle) = sin(angle) / cos(angle)

you note that e / c = sx * sin(alpha_x) / sx * cos(alpha_x)

The factor sx cancles leaving us with e / c = sin(aplha_x)/cos(alpha_x) = tan(alpha_x)

In the data example above e = 0 and c = 1.10401683181001. That means

e / c = tan(alpha_x) = 0. As a result we conclude that the x-axis rotation in this example is by an angle of 0 radians = 0 degree ==> no rotation at all.

Try the same for the y axis rotations. You are correct if you also find that there also was no y-axis rotation.

Since alpha_x and alpha_y are zero, c and f (multipliers of
cos(0)=1) show the respective scale changes from image to map along
the X and Y map coordinate axis = 1.1040168318 (final 4 digits
removed, we talk about this later).

This tells us that each pixel in the image transforms into map image pixel of about 1.1 m by 1.1 m in map coordinate extent.

The image concerned must have been an image previously rectified to
a UTM system (no rotations and scales close to 1). However, there
remains a slight difference, too large to be explained by map
projection difference (e.g. MGA on WGS94 aginst AMG on AGD66). The
scale differences then should not be larger than 0.0001 or 0.0002).
There must be another reason for these differences. Also remenber
that we finally dropped the last 4 digits of the scale values to
make them appear equal.

The reasons for this discrapancy are in the practical registration
process:

In order to solve an equation system with n unknown parameters (a to f, etc) you need exactly n independant equations to obtain a result. Meaning that we need to fill in X and Y map coordinates and x' and y' image coordinates for corresponding points (map and image) leaving just the parameters a to f unknown. Since each of these corresponding 'Match' or 'Control' Points deliver two of these equations (one for X = ... and one for Y = ...), we need to find at least two such points for the 4-parameter transformation or 3 for the 6-parameter transformation. (Note: These points should be selected as widely spaced as possible and in the 3 point case, not falling on a straight line, rather than forming the largest possible triangle).

The image registration process proceeds tehn as follows in most GIS systems:

Enter X and Y and x' and y' either manually or by mouse pointing.

ArcGIS allows pointing with the mouse for both, map and image points, while MapWindow required manual keyboard entry for map coordinates and pointing with the mouse for the corresponding image points.

In either case, the registration process creates **uncertainty**
in pointing to and clicking on the correct matching point(s).
Further more, when pointing to an image point, the exact pixel
address in x' and y' pixel counts is retrieved for the pixel
clicked on. There is only a very slim chance, that the
corresponding match point is also shown exactly at the center of a
pixel, which contributes a positional uncertainy of up 1/2 the
pixels size in both x' and y' (in case of the data example used
above, this could lead to about 0.5 m mis-registration per point
==> the apperent difference in scale somewhat points to 0.1 m
error)

Of course, we cannot exclude complete mis-interpretation of the
matching points in both, the map and the image (e.g. selecting the
wrong street intersection as control).

For this reason, it is good practice to allow the selection of more
than the required Control points, for example using 6 points
instead of 2 for the Similarity Transformation parameter
solution.

However, using just the algebraic approach in solving such an
'over-determined' equation system together with the inherent
uncertainties leads to multiple possible solutions rarely agreeing
with each other and leaving us with the uncertain decision of which
solution finally to use. Of course, we would like to obtain the
**most likely correct** parameter set with the **least
error** propagated from our uncertain input data.

The solution to this dilema is to apply statistically based
**Least Squares Adjustment** to the over-determined equation
system.

Least Squares adjustments allow (a) to detect extreme errors (e.g.
mis-pointing) and (b) a solution with the least errors for the
previously unknown parameters.

ArcGIS employs such Least Squares approach, at the time of writing this article MapWindow's Georeferencing plugin does not use Least Squares.

Finally, just for your exercise of interpreting parameter sets of
georeferenced images, here another example:

From a file by the name of *11-152-ortho_cut.jgw*

c = 0.983285005757208340

d = -0.00321150332591071470

e = -0.00445338734757580370

f = -0.986201839902504210

a = 480267.64268277911

b = 6978440.0628153477

Most GIS software packages follow a file name convention
(loosely).

The above example parameter file belongs to the image file
*11-152-ortho_cut.jpg*.

Note, the filenname before the '.' remains the same, while the letters after the '.' have slightly changed incorporating the letter 'w' for the parameter file. These parameter files are text files and can be viewed and edited using for example Notepad under Windows.

For more information search the web for 'Similarity Transformation', 'Affine Transformation', 'Least Squares Adjustment', Least Squares adjustment of observation equations', etc.

*This page may contain inconsistencies and errors. Please correct or complain to --Gngdowid 04:53, 25 December 2009 (UTC)*