Coordinate systems and reprojection

Every map file says where its cells are with coordinates, and every set of coordinates belongs to a coordinate reference system (CRS). Two maps of the same place in different systems do not line up; a map in the wrong system gives areas that are off by half. In this tutorial you will see the same landscape in degrees and in metres, measure a forest's area three ways and find that the naive way overestimates it by 60%, reproject a raster with raster.reproject, and learn why a land cover map must never be reprojected with the method that suits an elevation model.

The idea: from a round Earth to a flat map#

A position on Earth can be given in two ways.

• Geographic coordinates: longitude and latitude, in degrees, on a model of the Earth's shape (an ellipsoid and a datum that ties it to the ground; WGS 84, used by GPS, is the most common). They cover the whole globe, but a degree is not a fixed length: a degree of latitude is always about 111 km, while a degree of longitude shrinks from 111 km at the equator to nothing at the poles.
• Projected coordinates: x and y, usually in metres, on a flat plane. A map projection is the mathematical rule that flattens the curved surface onto the plane. On a projected map a cell of 30 m is 30 m wide everywhere, so areas and distances can be counted directly.

No projection is perfect: flattening a sphere always stretches something. A projection can keep areas (equal-area projections, such as Albers or Lambert azimuthal equal-area), or shapes and angles locally (conformal projections, such as Mercator and the Universal Transverse Mercator, UTM), or distances along some lines, but never all of them at once. Snyder's manual for the U.S. Geological Survey (Snyder, 1987) remains the classic description of the projections and their formulas; today they are computed by the PROJ library (PROJ contributors), which Terraflux uses through GDAL.

Names for coordinate systems: EPSG codes#

A CRS is described by a long text (a WKT), but most are known by a number in the EPSG registry, a public database of coordinate systems kept by the International Association of Oil & Gas Producers (IOGP). Terraflux accepts an EPSG code (EPSG:4326), a WKT or a PROJ string wherever it asks for a CRS. Some you will meet:

For a study area in one UTM zone, its UTM code is a good default for distances; for areas over larger regions use the continent's equal-area projection.

Why areas and distances need a projected CRS#

Here is the same question answered three ways: how much forest is on a forest map of 1000 x 1000 cells in degrees (landschaft's example habitat map; cells of 0.000305°, at about 51.4° N in central Germany), with 428,490 cells of forest?

The naive count treats every cell as a square, but at 51° N a cell of 0.000305° is only 21 m wide. Counting cells on an equal-area projection gives the right answer to 0.1%. The same goes for every measure that uses lengths or areas: patch areas, edge lengths, core areas, the distance thresholds of connectivity analyses, the dispersal distances of a species. Landscape metrics, MSPA's edge width and distances computed in degrees are wrong unless the method corrects for it.

When to reproject#

• To combine layers. Maps that are combined cell by cell must share one CRS (and one grid: see One grid for all your maps). Choose the CRS of your main map, usually the land cover, and bring the others onto it.
• To measure. Areas on an equal-area projection; local distances on UTM or another conformal projection.
• To share. Occurrence records, web maps and many global datasets use EPSG:4326.

Every reprojection of a raster resamples it: the new grid's cells do not fall on the old ones, so each new cell takes its value from the old cells near it. That always loses a little, so reproject each layer once, from its original file, rather than back and forth.

What you need#

The CRS of each of your files. Open Metadata in the viewer's header and read coordinate system: Terraflux shows its name and EPSG code (and the full WKT below). Files that say nothing (none) need you to tell Terraflux what they are; see A file without a CRS below.

Reproject a raster#

raster.reproject warps rasters to another CRS and/or resolution with GDAL's gdalwarp. Its settings:

• crs: the new CRS (EPSG:32610, a WKT, a PROJ string);
• resolution: the new cells' side, or [x, y], in the new CRS's units (30 for 30 m; about 0.0003 for degrees). Without it GDAL picks one that keeps about the same number of cells;
• resampling: how new cells get their values (next section);
• source_crs: the rasters' CRS, when their files do not say it;
• extent: [xmin, ymin, xmax, ymax] in the new CRS, to cut while warping;
• path: the file to write (else a file in a cache), and options for more of gdalwarp's options.

format = "terraflux/1"

[[nodes]]
id = "elevation"
op = "raster.read"
params = { path = "orelevation.tif" }
note = "SRTM elevation in EPSG:4326 (degrees), south-western Oregon."

[[nodes]]
id = "elevation_utm"
op = "raster.reproject"
inputs = { raster = "elevation.raster" }
params = { crs = "EPSG:32610", resolution = 30, path = "elevation_utm.tif" }
note = "UTM zone 10 N in metres; values interpolated (bilinear), the default for a continuous map."

[[plots]]
id = "plot"
title = "Elevation (UTM 10 N)"
layers = [{ from = "elevation_utm.raster" }]

The input is the raster output of raster.read, the file on disk: reprojection is done file to file, whatever the size of the raster, and the result is again a raster on disk.

Choose the resampling method by the data#

When a new cell falls between old cells, where does its value come from?

With resampling = "auto" (the default), Terraflux chooses by the data: for classes the nearest cell, or the most common one when the cells grow; for values bilinear, or the mean when the cells grow.

Why does it matter so much? Because interpolating class codes invents classes. Bilinear resampling of a cell between arable land (12) and coniferous forest (24) gives 18: pastures. Between discontinuous urban fabric (2) and arable land it gives 4 to 10: roads, ports, airports, dumps. Here is the 15 x 15 km window of the CORINE map around Hanover reprojected to EPSG:4326, once with near and once with bilinear:

near

bilinear

Over the whole window, bilinear gave 6.7% of the cells another class than nearest, and 4.5% of the cells a class that does not exist in the window at all: 26 invented classes, including rice fields, olive groves, glaciers and perpetual snow and salt marshes, in Lower Saxony. Landscape metrics computed on such a map would count thousands of tiny patches along every edge. With auto the result is identical to near.

Terrain in degrees#

Slope, aspect and hillshade compare heights in metres with horizontal distances. On an elevation model in degrees, the horizontal distances are in degrees too, unless the software corrects for it. Here is the mean slope of the Oregon elevation model computed by raster.terrain (GDAL's gdaldem) three ways:

The GDAL brought with Terraflux (3.12) recognises a surface in degrees and converts its distances by itself, so the default gives nearly the slope of the projected model. The option -s 111120 (one degree taken as 111,120 m in both directions), which older versions of GDAL needed, now underestimates the slopes at 43° N by about 10%, because it treats the shorter degrees of longitude as if they were degrees of latitude. The safest way, whatever the software, is to reproject the elevation model to a projected CRS in metres first.

A file without a CRS#

Some files have coordinates but do not say in which system. The land use of Sibuyan Island in the Philippines (see Your first raster) has x from 436,000 to 468,000 and y from 1,356,000 to 1,383,000, but no CRS in its header. Those numbers fit UTM zone 51 N (EPSG:32651), the zone of Sibuyan at 122.5° E. The setting source_crs says so:

format = "terraflux/1"

[[nodes]]
id = "landuse"
op = "raster.read"
params = { path = "sibuyan_landuse.tif" }
note = "No CRS in the file: its coordinates are those of UTM zone 51 N."

[[nodes]]
id = "wgs84"
op = "raster.reproject"
inputs = { raster = "landuse.raster" }
params = { crs = "EPSG:4326", source_crs = "EPSG:32651", path = "sibuyan_wgs84.tif" }
note = "Classes: the nearest cell (auto)."

The reprojected map falls at 122.41° to 122.71° E and 12.27° to 12.51° N: on Sibuyan. That check, the result landing where the place really is, is the way to confirm a guessed CRS; a basemap under the map (B in the viewer) shows it at a glance. raster.match has the same source_crs setting.

Common pitfalls#

• Layers that do not line up. Two maps of the same place that appear far apart, or one of them squeezed, are in different CRSs. Reproject one onto the other (better: onto the reference grid with raster.match).
• Degrees used as metres. A resolution of 30 given with crs = "EPSG:4326" asks for cells of 30 degrees. The resolution is in the new CRS's units.
• Latitude and longitude swapped. Occurrence records in a table have lat and lon columns; swapping them sends points to another continent, or into the sea.
• Measuring on Web Mercator. Areas inflated 2 to 4 times at mid-latitudes. Use an equal-area projection.
• A study area across UTM zones. UTM is accurate within its zone (6° of longitude); far beyond it, distortion grows. For large areas use an equal-area projection of the continent.
• Resampling classes with bilinear or cubic. Invented classes along every edge. Use near or mode.
• Different datums. WGS 84, ETRS89 and SIRGAS 2000 agree within about a metre; older datums (NAD27, Corrego Alegre, SAD69) can be tens of metres off. GDAL converts between them when the CRSs are known; it cannot when a file lacks its CRS.

Going further#

• One grid for all your maps: resolution, alignment, raster.match and the reference grid.
• Your first raster and Vector data.
• The operators: raster.reproject, raster.match, raster.terrain.

References#

• GDAL/OGR contributors. GDAL/OGR Geospatial Data Abstraction software Library. Open Source Geospatial Foundation. https://doi.org/10.5281/zenodo.5884351
• PROJ contributors. PROJ coordinate transformation software library. Open Source Geospatial Foundation. https://doi.org/10.5281/zenodo.5884394
• Snyder, J. P. (1987). Map Projections: A Working Manual (U.S. Geological Survey Professional Paper 1395). U.S. Government Printing Office. https://doi.org/10.3133/pp1395