Chapter 3 Polar maps

What is required for mapping in polar regions?

  • shape-feature orientations, the dateline and the poles
  • handling of data in long-lat, understanding of metrics and shape (area, length, angle)
  • projected maps can be very hard to do, each has its own issues
  • quantities like vector fields
  • conventions for global data -180/180, 0/360

3.1 Projections

A brief overview here. For more please see the geocompr book.

Projections are used to solve particular measurement problems. These include distance, area, and angle.

  • Equal-area - area is simple, we can calculate planar area anywhere
  • Equi-distant - less simple, applies along an axis or from one point only
  • Conformal - preserves shape, i.e. angles are sensible

No projection is perfect for every situation, and quite often tools will work in geographic coordinates, or raw longitude latitude values using algorithms suited for angular units on a sphere or ellipsoid.

Most projection specification is made by choosing an EPSG code, and if one exists for your region then it’s a good choice. However, sometimes the right EPSG does not exist and it’s simpler to specify a custom projection.

Each +proj= code is a projection family, there are many more than shown here.

Lambert Azimuthal Equal Area

+proj=laea +lon_0=0 +lat_0=-90 +datum=WGS84

Polar Stereographic (south)

+proj=stere +lon_0=0 +lat_0=-90 +lat_ts=-70 +datum=WGS84

Longitude-Latitude

+proj=longlat +datum=WGS84

Lambert Conformal Conic

+proj=lcc +lon_0=147 +lat_0=-42 +lat_1=-50 +lat_2=-35 +datum=WGS84

Each of these has in common a central longitude and latitude, and many custom projections can be made using only these parameters.

Other parameters include

  • +x_0, +y_0 - false easting (adds an arbitrary offset, usually to avoid negative coordinates)
  • +lat_0, +lat_2 - these are for conic projections, the secant latitudes
  • +ellps, +no_defs, +over - various datum/ellipsoid specific
  • +lat_ts - latitude of true scale projection specific, e.g. Stereographic
  • +zone, +south - projection-specific parameters, e.g. for UTM

Functions in R to expand EPSG parameters:

3.2 Mapping in R

## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /usr/share/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.3-1

The oldest general mapping tool in R is the maps package. It has a simple whole-world coastline data set for immediate use.

The data underlying this live map is available by capturing the output as an actual object.

If we look carefully at the southern edge and the eastern edge, notice that the coastline for Antarctica does not extend to the south pole, and the Chukotka region of Russia east of 180 longitude is not in the western part of the map.

A very similar and slightly more modern data set is available in the maptools package.

(It’s not sensible use maps or maptools for coasline or country boundary data, but they are handy for exploring concepts ).

This data set aligns exactly to the conventional -180/180 -90/90 extent of the longitude/latitude projection.

3.3 Exercise 1

How can we find the longitude and latitude ranges of the maps data maps_c and the maptools data wrld_simpl?

## List of 4
##  $ x    : num [1:82403] -69.9 -69.9 -69.9 -70 -70.1 ...
##  $ y    : num [1:82403] 12.5 12.4 12.4 12.5 12.5 ...
##  $ range: num [1:4] -180 190.3 -85.2 83.6
##  $ names: chr [1:1627] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda" ...
##  - attr(*, "class")= chr "map"

3.3.1 EX 1 ANSWER

SOLUTION 

## [1] -180.0000  190.2708
## [1] -85.19218  83.59961
## [1] -180.00000  190.27084  -85.19218   83.59961

3.4 Exercise 2

Can we draw polygons with a fill colour with the maps package?

Why don’t these work?

3.4.1 EX 2 ANSWER

SOLUTION

The actual data structure in the maps package is not suitable generally for subsetting polygons, or controlling polygon drawing.

3.5 Polygon catastrophe

What’s going on? Look at the very south-eastern corner of the map. The “coastline” has been extended to the very south boundary of the available area.

When we add the old maps coastline see that it does not extend to 90S and it does not traverse the southern boundary.

One reason for this is that if we choose a projection where the east and west edges of the Antarctic coastline meet then we get what looks a fairly clean join.

If we try the same with wrld_simpl it’s not as neat. We have a strange “seam” that points exactly to the south pole (our projection is centred on longitude = 0, and latitude = -90.

3.6 Exercise 3

How can we add our own data, a set of points, onto the polar map?

Add the pts matrix of longitude/latitude to this map, or define your points and add those.

Hint: we need the coordinate reference system string “+proj=laea +lat_0=-90 +datum=WGS84” and a function to transform a matrix of longitude-latitude values.

3.7 Let’s use the maps data!

In maps_c we have the maps data structure, and this looks promising.

## List of 4
##  $ x    : num [1:82403] -69.9 -69.9 -69.9 -70 -70.1 ...
##  $ y    : num [1:82403] 12.5 12.4 12.4 12.5 12.5 ...
##  $ range: num [1:4] -180 190.3 -85.2 83.6
##  $ names: chr [1:1627] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda" ...
##  - attr(*, "class")= chr "map"
## [1] -12709814  12704237 -12576156  12470787

The problem is that the maps database has enough internal structure to join lines correctly, with NA gaps between different connected linestrings, but not enough to draw these things as polygons. A similar problem occurs in the default projection. While wrld_simpl has been extend by placing two dummy coordinates at the east and west versions of the south pole, this data set does not have those.

We have to look quite carefully to understand what is happening, but this is wrapping around overlapping itself and so close to the southern bound we barely notice.

## [1] -20037508  20037508 -20179524  18351859

3.8 Reprojecting to polar regions

In essence it should be easy, but details really matter. Here the details of the format, the orientation and layout of the data (NetCDF, and common conventions used in physical ocean models)

## Loading required namespace: ncdf4
## Warning in .varName(nc, varname, warn = warn): varname used is: sst
## If that is not correct, you can set it to one of: sst, anom, err, ice
## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
## 48 projected point(s) not finite
## Warning in rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], :
## 6211 projected point(s) not finite

## Warning in rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, :
## 48 projected point(s) not finite
## Warning in rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], :
## 6213 projected point(s) not finite

Straight lines in longitude/latitude do not maintain their shape when reprojected unless we add extra vertices along lines of constant longitude and latitude. The longitudes are great circles, but the latitude lines are not.

The pole is undefined in Mercator projection.

##           [,1]         [,2]
## [1,] -20037508 7.081155e-10
##           [,1]     [,2]
## [1,] -20037508 -8362699
## Warning in rgdal::project(cbind(-180, -90), "+proj=merc"): 1 projected
## point(s) not finite
##      [,1] [,2]
## [1,]  Inf  Inf

Sensible polar projections are stereographic and Lambert Azimuthal Equal Area. These are more or less identical over a large area, so knowing what is in use really matters!

This video explains the Lambert Azimuthal Equal Area projection: https://www.youtube.com/watch?v=quzIU4nL9ig (crux at 2:25 with the spherical shells).

3.9 Demonstration of the projection concept

3D plots that show various projection concepts.

First prepared here: http://mdsumner.github.io/2016/01/26/Three_Projections.html