Coordinates

This is a back-up of the WIKI.
We're working on a new wiki.

Main Page | Recent changes | View source | Page history | Log in / create account |

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:

1. Polar coordinates: involving one or two angles and a radial distance
2. Rectangular cartesian coordinates

Polar Coordinates [1]

Points in the (2D) polar coordinate system with pole O and polar axis L. In green, the point with radial coordinate 3 and angular coordinate 60 degrees, or (3,60??). In blue, the point (4,210??).

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

2D: and 3D:

For mapping and GIS applications the following is preferred:

A. 2D coordinates (x,y)

1. Plane x, y coordinates (dominant case):
• x and y axis are straight ordinate axis spanning open a plane surface
2. (curved) surface x, y coordinates:
• x and y axis are straightest ordinate axis on a curved surface (ie. on a globe)
3. x and y axis intersect in a point called the 'Origin' of the coordinate system
4. x and y axis intersect at right angle
5. distances along each axis are referenced against this origin : x=0 and y=0
6. distance scale is set either as linear (metric; preferred case) or as a non-linear function of the distance itself (e.g. logarithmic)
7. 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)
8. (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
• 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
• 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.
• 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:

1. Operating Environment: GUI (Graphical User Interface), Data I/O (basic file access and screen display) Program language support (libraries)
2. Database Management System (DBMS): SQL, MYSQL, Postgre, etc. for programmed and user defined data access and interogation
3. GIS specific Transformations of positional data
4. 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.