Projection support in QGIS

We are in the process of implementing support for Projections in QGIS. This document provides a rough guideline to the approach we are taking to achieve this.

Please note this document is still under construction

Tim Sutton Jan 2005

Contents

  1. Projection support in QGIS
    1. Key Concepts
    2. Behaviour Rules
    3. Changes to QGIS Source Required
      1. SRID
      2. Calculating Project Extents
      3. The render process
    4. The Projection Databases
    5. Comparing Projections
    6. IRC Discussions that need to be tidied up and included in this doc
    7. TODO:
    8. Useful bookmarks & further reading

Key Concepts

  1. On-the-fly: QGIS will be able to load multiple layers from different coordinate systems and reproject them into a common project coordinate system as the layer is rendered. This means that it is not neccessary to pre-process layers and save them to disk in the common coordinate system.

  2. Project Coordinate System (PCS): This is the common coordinate system that all layers will be reprojected to (if needed). Having a common coordinate system allows us to overlay multiple files in different coordinate systems into a single map display. In the source code, the project coordinate system is often referred to as the Dest (destination) coordinate system.

  3. Layer Coordinate System (LCS): Projection in QGIS will be layer centric. This means that each layer will know what its own Coordinate System (CS) is, and what the project level CS is. In the source code, the layer coordinate system is often referred to as the Source coordinate system.

  4. Well Known Text (WKT): This is a standard format (described in detail here) promoted by the OpenGis consortium for describing a coordinate system.

  5. SRID: This is an alternative approach for sepcifying a CS, used by the Postgis folks and described in detail in the OpenGIS Consortium Simple Features specification

  6. Spatial Reference System (SRS): A projected or unprojected spatial referencing system.

  7. FORWARD Transform: A CS transformation from the layer's native CS to the project CS

  8. INVERSE Transform: A CS transform from the common project CS to the layers CS

Behaviour Rules

http://mrcc.com/~gsherman/proj_set_layer_cs.png

Changes to QGIS Source Required

We are trying to make the Projection implementation as minimally invasive as possible. The logic relating to performing a

FORWARD or INVERSE transformation is encapsulated within the QgsCoordinateTransform class. Because of the prevalence of the WKT standard, the primary means for for defining the Source and Dest CS is by neabs if the QgsCoordinateTransform constructor. After the QgsCoordinateTransform object has been instantiated, the project coordinate system can be redefined using 

void setDestWKT(QString theWKT);

QgsMapLayer (which is the superclass for all layer types) now has a member to hold the QgsCoordinatSystem: 

QgsCoordinateTransform * mCoordinateTransform;

In addition, QgsMapLayer now defines the following virtual method: 

virtual QString getProjectionWKT()  = 0 ;

Which must be reimplemented by any class derived from the QgsMapLayer superclass. The provider architecture provides a specialised implementation class for each data format supported e.g. OGR, GPS, GRASS, POSTGIS etc. Currently the raster implementation in QGIS does not have a provider like archtecture, and so determining the WKT of a raster layer is handled directly in the QgsRasterLayer class. In the case of vector layers, The QgsVectorLayer class however does use this provider implementation, and thus determining the WKT of a layer needs to be delegated down to the provider instance associated with that layer: 

QString QgsVectorLayer::getProjectionWKT()
{
  //delegate to the provider
  if (valid)
  {
    return dataProvider->getProjectionWKT();
  }
  else
  {
    return NULL;
  }
}

Each provider needs to 'know' how to obtain the WKT for the datasets it can read, In some cases such as the GPS provider, the getProjectionWKT is hardcoded to return GEOCS/WGS84 (the gps babel software used internally by the GPS provider specifies that all gps layers are in this CS). In other cases (such as reading from ESRI .prj files). The WKT datum returned may not always be readable. For this reason a simple check is made during the QgsCoordinateSystem initialisation for the presence of a DATUM node, and if its not available the OGR myInputSpatialRefSys.morphFromESRI(); functions is used to 'clean up' the datum.

For the most part we are using the OGR Spatial Reference management tools to deal with the projection definitions. We tried using the OGR Coordinate Transformation tool to perform the actual reprojection operation but ran into problems when the reprojection involved a datum shift (for example from NAD27 to NAD83). For this reason the actual projection logic (at the vertex level) is being implemented directly using the proj library (OGR Coordinate Transformation is a wrapper for proj anyway!). Because we are using proj, setting the source and destination WKT of a QgsCoordinateTransform object need to be converted to proj format whenever they are changed and stored in the QgsCoordinateTransform private members mDestinationProjection and mSourceProjection as shown in the code snippet below. 

  char *proj4dest;
  myOutputSpatialRefSys.exportToProj4(&proj4dest);
  // store the dest proj parms in a QString because the pointer populated by exportToProj4
  // seems to get corrupted prior to its use in the transform
  mProj4DestParms = proj4dest;
  // check to see if we have datum information -- if not it might be an ESRI format WKT
  if(mProj4DestParms.contains("datum") == 0)
  {
    // try getting it as an esri format wkt
    myOutputSpatialRefSys.morphFromESRI();
    myOutputSpatialRefSys.exportToProj4(&proj4dest);
    mProj4DestParms = proj4dest;
  }
  if(mProj4DestParms.contains("datum") == 0)
  {
    throw QgsCsException("Destination coordinate system does not contain datum. WKT is : " + mDestWKT );
  }
  // init the projections (destination and source)
   mDestinationProjection = pj_init_plus((const char *)mProj4DestParms);
   mSourceProjection = pj_init_plus((const char *)mProj4SrcParms);

These members are pre-initialised proj transformers so they are readily available whenever a vertex needs to be transformed.

SRID

Spatial Reference System Identifiers (SRID's) are another popular format for uniquely expressing a Spatial Reference System (SRS) definition. The identifier is a unique serial number given to each SRS. In general this number will be identical to the EPSG number for the SRS. Confused yet? Yes there are a number of competing formats for describing SRS's. SRID's are of interest to us because they provide access to a large (~1800) collection of predefined projection definitions. The information associated with an SRID includes the WKT definition too so they are not incompatible with the implementation plans described above.

Having access to such a large colelction of defined WKTs is an opprotunity we do not want to miss out on - especially as the SRID database is independently maintained so its relatively easy for us to provide a current, accurate and up to date list of coordinate systems with each release of QGIS. This database will be used whenever the user is asked to make a selection of a coordinate system in a dialog in the QGIS gui. However there is one snafu....the OGR tools do not provide a ready way in which to determine the SRID given a WKT definition. This functionality is needed so that we can highlight the current projection of a layer in the aforementioned dialog boxes. A work around to this problem is to use the isSame() OGR function to iterate through the SRID database, retrieve its WKT and compare it to the layer WKT or the project WKT as appropriate. This may involve a fair bit of processing so we need to investigate if there is a more optimal way to achieve this.

Calculating Project Extents

Whenever a new map layer is added to the project, the canvas needs to recalculate the full extent of all layers in the project. The void QgsMapCanvas::recalculateExtents() (API) method iterates through each layer in the map layer registry and asks it for its extent. The extent is returned in the layer's native projection and reprojected into the project project so that it can be compared with the extents of all layers in the project.

The render process

http://mrcc.com/~gsherman/proj_draw_layers.png

When a user zooms or pans the map display, or changes the map extents in some other way (e.g. zoom to all), the current visible area extent for each layer need to be calculated in the native coordinate system of the layer. In order to do this an inverse transform needs to be carried out which will convert the mapcanvas extents (stored in the project coordinate system) into the layer coordinate system. This happens within the main render loop of the map canvas: 

void QgsMapCanvas::render(QPaintDevice * theQPaintDevice)

Within the main render loop, the registry is once again iterated through and every layer that is visible (i.e has rendering enabled) is passed the inverse projected project extents to its draw() method. This inverse projected extent is then used to select just the features what are within the map view to be drawn. If the coordinate systems of the layer and project match, the inverse transform is 'short circuited' and the project extents are simply returned unchanged.

Once the area of interest has been determined for the current layer being rendered and the appropriate selection of features has been made, the draw method starts iterating through each feature (in the case of a vector layer). As feature is retrieved, it is 'dissolved' down to its constituant vertices. Each vertex needs to be converted from map coordinate space to pixel coordinate space in order for it to be placed correctly on the screen. This is a three step process:

  1. Retrieve coordinates of vertex in layer's native coordinate system
  2. Forward transform the vertex into the common project coordinate system
  3. Convert the vertex from project coordinate system to map canvas pixel coordinates

The Projection Databases

After some conseideration we elected to use the sqlite3 file based database for persistant storage of projection definitions. SQLite is a lightweight, supported on all major operating systems and very efficient. This efficiency was a major consideration as accessing projection catalogues stored in simple text files proved to be a major overhead. SQLite databases are given the extension .db.

There are three databases associated with projection catalogue management:

The system spatial reference system dabatase

${PREFIX}/share/qgis/resources/srs.db (where ${PREFIX} denotes the QGIS installation prefix)

The users harvested SRS

~/.qgis/srs.db

The users custom projection definitions

~/.qgis/customproj.db

Comparing Projections

We are electing to use proj parameter strings as the lowest common denominator for internal (to QGIS) representation of projections. In order to compare two projections, we will follow the following procedure:

proj_workflow_equivalence_test

IRC Discussions that need to be tidied up and included in this doc

The coordinate system is represented by class QgsSpatialRefSys

It is the job of this class to store the reference system and to instantiate a ref sys using one of the folloowing standards:

Internally projection is done using libproj4 and so when QgsSpatialRefSys is initialised, the goal is to populate mParameters member with a valid proj4 string The actual coordinate transform is carried out by class QgsCoordinateTransform (xfrm for short from now on)

xfrm has two QgsSpatialRefSys members

I say 'usually' because of course you could use a xfrm object to carry out arbitary cs translations not related to layers At the moment each provider either honours one of the following calls:

(method names from memory, may be slightly different in real life)

One of the things I am busy doing is rearchitecting the proj system so that except for the data facing (i.e. provider classes) qgis is only aware of the QgsSpatialRefSys concept and not of SRID or WKT as currently implemented above

TODO:

Useful bookmarks & further reading

Projection support in QGIS (last edited 2005-05-09 00:25:40 by MateuszLoskot)