Chapter 3 Projections and geodesic calculations
If only every study fit within a state plane zone or UTM zone… but alas, they do not.
How can we check the coordinate reference system and transform coordinates into another system? How can we calculate area and length geodesically, rather than planimetrically?
3.1 Projections
Let’s check the projection of the watersheds layer.
## [1] "NAD83(2011) / Rhode Island"And let’s see what the geometry data looks like.
## [1] "MULTIPOLYGON(((107710.4025421 90042.7079543,108460.4025421 89292.7079543,108460.4025421 88292.7079543,107710.4025421 88292.7079543,107710.4025421 90042.7079543)))"
## [2] "MULTIPOLYGON(((108460.4025421 88792.7079543,108960.4025421 88792.7079543,109210.4025421 89042.7079543,109460.4025421 89042.7079543,109460.4025421 88292.7079543,108460.4025421 88292.7079543,108460.4025421 88792.7079543)))"
## [3] "MULTIPOLYGON(((109460.4025421 88292.7079543,109460.4025421 89042.7079543,109210.4025421 89042.7079543,108960.4025421 88792.7079543,108460.4025421 88792.7079543,108460.4025421 89292.7079543,107710.4025421 90042.7079543,107710.4025421 90292.7079543,109710.4025421 90292.7079543,109710.4025421 88292.7079543,109460.4025421 88292.7079543)))"Now let’s transform the data into NAD 1983 geographic coordinates.
And check the CRS.
## Coordinate Reference System:
## User input: EPSG:4269
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Geodesy."],
## AREA["North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands."],
## BBOX[14.92,167.65,86.45,-40.73]],
## ID["EPSG",4269]]And inspect the geometry data to see if the coordinates have changed.
## [1] "MULTIPOLYGON(((-71.4070897 41.8940335,-71.3980629 41.8872734,-71.3980772 41.8782701,-71.4071125 41.8782778,-71.4070897 41.8940335)))"
## [2] "MULTIPOLYGON(((-71.3980701 41.8827718,-71.3920461 41.8827663,-71.3890303 41.8850142,-71.3860182 41.8850113,-71.3860302 41.8782588,-71.3980772 41.8782701,-71.3980701 41.8827718)))"
## [3] "MULTIPOLYGON(((-71.3860302 41.8782588,-71.3860182 41.8850113,-71.3890303 41.8850142,-71.3920461 41.8827663,-71.3980701 41.8827718,-71.3980629 41.8872734,-71.4070897 41.8940335,-71.4070864 41.8962844,-71.3829856 41.8962623,-71.3830185 41.8782558,-71.3860302 41.8782588)))"Indeed, they are now stored as latitude and longitude coordinates!
3.2 Geodesic calculations
As it turns out, sf will calculate spherical lengths and areas with st_area and st_length if the geometries are stored in a geographic coordinate system.
This is tested with the st_is_longlat() function.
Furthermore, there is a global setting for sf to switch spherical calculations off and use more accurate ellipsoidal calculations. This can be achieved with sf_use_s2(FALSE).
For perimeter calculations below, st_perimeter(geom) should work, and does work when running code blocks manually.
For reasons I cannot explain, the operation needs to be simplified to a length function on a line in order for me to compile this code into a book.
sf_use_s2(TRUE)
watersheds_NAD83 <- watersheds_NAD83 |> mutate(
spherical_area = st_area(geom),
spherical_perim = st_length(st_cast(st_cast(geom, "POLYGON"), "LINESTRING"))
)
sf_use_s2(FALSE)
watersheds_NAD83 <- watersheds_NAD83 |> mutate(
ellipsoidal_area = st_area(geom),
ellipsoidal_perim = st_length(st_cast(st_cast(geom, "POLYGON"), "LINESTRING"))
)
watersheds_NAD83 |> st_drop_geometry() |> kable()| name | spherical_area | spherical_perim | ellipsoidal_area | ellipsoidal_perim |
|---|---|---|---|---|
| Moshassuck | 1029719.6 [m^2] | 4561.014 [m] | 1031261.2 [m^2] | 4560.685 [m] |
| Narragansett | 592869.2 [m^2] | 3350.138 [m] | 593756.2 [m^2] | 3353.571 [m] |
| Blackstone | 2371471.7 [m^2] | 7909.285 [m] | 2375025.2 [m^2] | 7914.255 [m] |
Finally, you cannot calculate geodestic distance and area with geometries using projected coordinate systems in sf.
You could achieve this by transforming temporarily as you calculate, e.g. st_area(st_transform(geom, 4269)).
For example, we can calculate the planar and ellipsoidal areas of watersheds as follows:
sf_use_s2(FALSE)
watersheds_measure <- watersheds |>
mutate(planar_area = st_area(geom),
ellipsoidal_area = st_area(st_transform(geom, 4269)))
watersheds_measure |> st_drop_geometry() |> kable()| name | planar_area | ellipsoidal_area |
|---|---|---|
| Moshassuck | 1031250 [m^2] | 1031261.2 [m^2] |
| Narragansett | 593750 [m^2] | 593756.2 [m^2] |
| Blackstone | 2375000 [m^2] | 2375025.2 [m^2] |