The east coast of the United States spans nearly 20 degrees of longitude and 30 degrees of latitude, with a major axis that runs more “northeast-to-southwest” than “north-to-south”. Because of this, I wind up with a lot of wasted space when trying to plot the coast. Take a look to see what I mean:
library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
library(sf)
Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(rnaturalearth)
Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(ggplot2)theme_set(theme_minimal())
east_coast <-# Import the Natural Earth states polygon using rnaturalearthne_states(country ='united states of america',returnclass ='sf') |># crop to the US east coastst_crop(xmin =-82, xmax =-67,ymin =20, ymax =48)
The rnaturalearthhires package needs to be installed.
── R CMD build ─────────────────────────────────────────────────────────────────
* checking for file 'C:\Users\darpa2\AppData\Local\Temp\RtmpmmpUkt\remotesc40603f1bc5\ropensci-rnaturalearthhires-db8e433/DESCRIPTION' ... OK
* preparing 'rnaturalearthhires':
* checking DESCRIPTION meta-information ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building 'rnaturalearthhires_1.0.0.tar.gz'
Installing package into 'C:/Users/darpa2/AppData/Local/R/win-library/4.3'
(as 'lib' is unspecified)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
ggplot() +geom_sf(data = east_coast)
Firstly, I maybe should choose a different spatial polygon – this one has sunk the entirety of Virginia’s eastern shore! I’ll ignore that for now and move on.
If I’m planning to map migrations of fish along the coast (a common occurrence), I only need a little bit along the coast. The Sargasso Sea to the bottom right and Great Lakes to the top left are completely unnecessary.
One way to address this is to rotate the map so that the axis running from Florida to Maine is more-or-less vertical. My first attempt at this was using an affine transformation as outlined in the sf package documentation, but that didn’t quite work out.
an affine transformation…is a geometric transformation that preserves lines and parallelism, but not necessarily Euclidean distances and angles.
In addition to the sf vignette linked above, the free Geocomputation with R book provides a good outline. I’ll use their rotating function below to rotate the polygon counterclockwise by 40 degrees, which should make that N-S axis.
# Rotation functionrotation <-function(a){ r <- a * pi /180# degrees to radiansmatrix(c(cos(r), sin(r), -sin(r), cos(r)), nrow =2, ncol =2)} # Apply affine transformation by rotating the geometry columnec_affine <- east_coast |>mutate(geometry = geometry *rotation(-40))ggplot() +geom_sf(data = ec_affine)
Alright! Now all we need to do is trim the plotting area down and we’re in business!
But… wait a minute. Those axes make no sense. Since when is Maine at -10 degrees latitude? Let’s compare how the coordinates were changed.
It looks like the raw values were changed (as anticipated). However, they’re values that don’t really make sense in terms of the coordinate reference system with which we started. There is also no longer a CRS listed! So, our polygon is “spatial” insofar as it having a spatial attributes via a well-known text geometry column as made by sf, but it’s no longer actually something that we know how to place in space.
This is where the problems begin.
What if, for instance, we wanted to add a layer from a different source?
Error in st_transform.sfc(X[[i]], ...): cannot transform sfc object with missing crs
We get a blank plot and an error – since the affine-transformed polygon no longer has a CRS, sf and ggplot no longer know where or how to place the new object onto the map.
The difference between the Mercator (what we see when we look at something like Google Maps…ish; that’s actually web Mercator), transverse Mercator (the “TM” of “UTM” fame), and oblique Mercator lies in what it considers to be the major axis. A Mercator uses the equator, while a transverse Mercator uses a meridian (a particular longitude). Oblique Mercator uses some arbitrary line. This line can be determined using a point (the +lonc and +lat_0 arguments in a PROJ string) and an azimuth (either the +alpha or +gamma) or two points (+lat_1, +lon_1, +lat_2, +lon_2). I’ve found that providing a rotation away from North in degrees to the +gamma argument is what makes the most-intuitive sense to me.
This seems to be doing what we want at first blush. The graticule is rotated in a way that makes sense, and the axes are consistent with what we would think the proper values would be. What does it look like under the hood?
Well, that…didn’t work. Why not? Well, in essence, we’re providing two different units: ec_omerc has units of meters and we’ve provided units of degrees to coord_sf, assuming that it’s still using latitude and longitude (a CRS of WGS83/EPSG4326 in this case).
The hint to solve this is in the help documentation of coord_sf, specifically in the default_crs argument:
The default CRS to be used for non-sf layers (which don’t carry any CRS information) and scale limits. The default value of NULL means that the setting for crs is used. This implies that all non-sf layers and scale limits are assumed to be specified in projected coordinates. A useful alternative setting is default_crs = sf::st_crs(4326), which means x and y positions are interpreted as longitude and latitude, respectively, in the World Geodetic System 1984 (WGS84).
When we use the xlim and ylim arguments, it’s going to default to using the CRS that is in the first layer. This would be our oblique Mercator projection, so we would have to provide something in meters. It’s so easy to just look at the axes that have been printed and work from there – so, we’re thinking in WGS84/EPSG 4326. We can just tell coord_sf that this is our frame of reference using the default_crs argument.
Almost there. But, wait – shouldn’t North Carolina be in the picture if we’re down around 37 degrees latitude? Let’s check using the same limits in the original projection.
Yes – North Carolina is there. So what’s going on?
Well, we are seeing 37 degrees latitude in the oblique Mercator projection – it’s the little sliver on the bottom right. If we’re rotating everything counter-clockwise by 40 degrees and keeping the same rectangular view, this means that we’re adding some latitude at the bottom right of the frame and removing some from the top left. Similarly, we’d add some longitude in the bottom left and remove some from the top right. We need to fudge things now to get what we want included.
If you want to rotate your map, you can do it in at least two ways:
Affine transformation
Oblique Mercator
If you use the Affine transformation, you’re just multiplying you coordinates by a number. This will rotate everything and make it plot-able, but you’re now using a custom, undefined coordinate reference system.
If you project the underlying spatial data into some sort of rotated coordinate reference system, like the oblique Mercator outlined here, you have a defined coordinate reference system and can continue to treat your polygon as a spatial object
When plotting with ggplot2, you can continue to use intuitive latitude/longitude to “zoom in”, but you will likely have to provide a default CRS when the units are not in degrees
Using a rotated projection and a rectangular view port will likely produce unintuitive results – be prepared to work iteratively and fudge things around.
