Skip to main content

R-GIS for Ecology

2 Coordinate Systems

2.1 Learning objectives

  1. Understand the basic concepts of Coordinate Reference Systems (CRS)

  2. Distinguish between geodetic (geographic) and projected CRS

  3. Learn how to select an appropriate CRS based on data type and analysis goals

In this chapter, we will work with the following packages. Before starting any exercises, make sure to run the following code.

if (!require(pacman)) install.packages("pacman")

pacman::p_load(tidyverse,
               sf,
               mapview)

2.2 Getting started

What are Coordinate Reference Systems (CRS)? CRS define how spatial data are projected onto the Earth’s surface—essentially, how locations on the curved Earth are represented on a flat map.

To get a sense of how CRS works in practice, let’s explore a dataset on fish distributions from North Carolina. This regional dataset contains georeferenced sampling sites using longitude and latitude coordinates.

We’ll start by reading the dataset from the data sub-directory in your R project:

df_fish <- read_csv("data/data_finsync_nc.csv")

This dataset contains the following columns:

Column Description
site_id Site name
year Sampling year
date Sampling date
lon Longitude
lat Latitude
latin Scientific name

Let’s take a quick view of the dataframe:

print(df_fish)
## # A tibble: 1,426 × 6
##    site_id               year date         lon   lat latin                  
##    <chr>                <dbl> <date>     <dbl> <dbl> <chr>                  
##  1 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Clinostomus funduloides
##  2 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Etheostoma flabellare  
##  3 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Hybopsis hypsinotus    
##  4 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Lepomis auritus        
##  5 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Moxostoma rupiscartes  
##  6 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Nocomis leptocephalus  
##  7 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Notropis chiliticus    
##  8 finsync_nrs_nc-10013  2009 2009-06-27 -81.5  36.1 Semotilus atromaculatus
##  9 finsync_nrs_nc-10014  2009 2009-06-03 -80.4  35.9 Clinostomus funduloides
## 10 finsync_nrs_nc-10014  2009 2009-06-03 -80.4  35.9 Fundulus rathbuni      
## # ℹ 1,416 more rows

Since this dataframe contains longitude and latitude columns, we can make this data as georeferenced. To do so, we’ll need to use functions from sf package, which makes the mere numbers of coordinates linked to CRS of your choice. Specifically, we first convert the original dataframe to unique combinations of logitude & latitude (= unique location) with dplr::distinct() then use sf::st_as_sf():

sf_site <- df_fish %>% 
  distinct(site_id, lon, lat) %>% # get unique combinations of longitude & latitude
  st_as_sf(coords = c("lon", "lat"),
           crs = 4326)

print(sf_site)
## Simple feature collection with 122 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.02452 ymin: 34.53497 xmax: -76.74611 ymax: 36.54083
## Geodetic CRS:  WGS 84
## # A tibble: 122 × 2
##    site_id                          geometry
##  * <chr>                         <POINT [°]>
##  1 finsync_nrs_nc-10013 (-81.51025 36.11188)
##  2 finsync_nrs_nc-10014 (-80.35989 35.87616)
##  3 finsync_nrs_nc-10020 (-81.74472 35.64379)
##  4 finsync_nrs_nc-10023  (-82.77898 35.6822)
##  5 finsync_nrs_nc-10024 (-77.75384 35.38553)
##  6 finsync_nrs_nc-10027 (-83.69678 35.02467)
##  7 finsync_nrs_nc-10029 (-80.65668 35.16119)
##  8 finsync_nrs_nc-10034 (-82.04497 36.09917)
##  9 finsync_nrs_nc-10041 (-80.46558 35.04365)
## 10 finsync_nrs_nc-10049 (-77.96322 34.58249)
## # ℹ 112 more rows

It looks similar to what we had as a tibble, but this version includes additional lines of metadata that describe key geospatial information. This data frame now represents a point layer (see the Geometry type), and includes a bounding box that defines the longitudinal and latitudinal extent (Bounding box). It also specifies the CRS in use (Geodetic CRS). The object class now includes sf.

class(sf_site)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"

Once the data frame is converted to the sf class, you can apply a variety of spatial mapping functions available in R. For example, we can quickly visualize the data using the mapview::mapview() function:

mapview(sf_site,
        legend = FALSE)

Fantastic — it works! This marks an important first step toward harnessing the full potential of RGIS.

You may save this object using saveRDS() function in R – in doing so, you can load this object for later use (see Chapter 3):

saveRDS(sf_site,
        file = "data/sf_finsync_nc.rds")

But, heads-up! If you don’t fully understand the CRS in use, there’s a high risk of performing misinformed analyses — potentially serious enough to invalidate your conclusions. Therefore, it’s important to take the time to understand several key aspects of CRS before conducting any spatial analysis. One of those is the difference between Geodetic vs. Projected CRS. The CRS we used in the fish data example is WGS 84, which is a type of geodetic CRS.

2.3 Geodetic vs. Projected

CRS are essential for accurately representing spatial data on the Earth’s surface. A geodetic CRS models the Earth as a 3D ellipsoid and represents locations using angular coordinates — latitude and longitude. These systems are ideal for global positioning and mapping where preserving true geographic coordinates is important. However, because they represent curved surfaces on a spherical model, they are not suitable for accurate distance or area calculations on flat maps. Common examples of geodetic CRS include WGS 84 (EPSG:4326), which is used by GPS and most global datasets, and NAD83 (EPSG:4269), which is widely used in North America.

In contrast, a projected CRS transforms geographic coordinates onto a 2D surface using mathematical projections. This allows for accurate measurements of distance, area, or direction in localized regions. Projected systems distort some properties (e.g., shape, area, or scale) but can be tailored for specific mapping needs. Examples of projected CRS include UTM zones (e.g., UTM Zone 17N, EPSG:32617) and State Plane systems like North Carolina State Plane (EPSG:32119), which are optimized for regional accuracy and are commonly used in government or environmental datasets.

Visit this website for visual aids.

Comparison of geodetic and projected CRS.
CRS Type Representation Advantages Disadvantages
Geodetic Latitude and longitude on an ellipsoidal model of Earth

Global coverage (e.g., GPS)

Standard for data exchange and web maps

Easy to interpret (lat/lon)

NOT for distance, area, or direction calculations

Curved surface not ideal for flat map operations
Projected X/Y Cartesian coordinates on a flat surface derived from a projection of the Earth

Enables accurate measurements (distance, area)

Optimized for regional analysis

Supports specialized map projections

Distorts some geographic properties

Limited to specific regions

Must choose the right projection for your purpose

2.4 CRS transformation

What if we want to perform spatial analysis on data provided in a geodetic CRS, like in the fish data example? Can we still calculate distances or areas? At first glance, it might seem impossible — but no worries. Thanks to extensive work in geospatial science, we now have well-established methods for transforming between CRS.

For example, suppose we want to calculate the distance between two sampling sites — say, the first two rows in our data frame.

sf_ft_wgs <- sf_site %>% 
  slice(c(1, 2))

print(sf_ft_wgs)
## Simple feature collection with 2 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -81.51025 ymin: 35.87616 xmax: -80.35989 ymax: 36.11188
## Geodetic CRS:  WGS 84
## # A tibble: 2 × 2
##   site_id                          geometry
##   <chr>                         <POINT [°]>
## 1 finsync_nrs_nc-10013 (-81.51025 36.11188)
## 2 finsync_nrs_nc-10014 (-80.35989 35.87616)

Let’s walk through how to do this using a projected CRS. If we want to calculate the distance between these two sites, we can use the sf::st_distance() function — but only after transforming the coordinates to a projected CRS1. Here, we’ll use UTM Zone 17N, which is well-suited for projecting spatial data in North Carolina. To do this, we’ll apply sf::st_transform() to convert our geodetic CRS (e.g., WGS 84) into the projected CRS for UTM Zone 17N based on WGS 84 (EPSG:32617).

sf_ft_utm <- sf_ft_wgs %>% 
  st_transform(crs = 32617)

print(sf_ft_utm)
## Simple feature collection with 2 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 454077.2 ymin: 3970402 xmax: 557782.3 ymax: 3996478
## Projected CRS: WGS 84 / UTM zone 17N
## # A tibble: 2 × 2
##   site_id                        geometry
## * <chr>                       <POINT [m]>
## 1 finsync_nrs_nc-10013 (454077.2 3996478)
## 2 finsync_nrs_nc-10014 (557782.3 3970402)

Now, the CRS metadata shows: Projected CRS: WGS 84 / UTM zone 17N. You’ll also notice that the bounding box values have changed to large numbers (check with st_bbox(sf_ft_utm)), which clearly are not in longitude or latitude.

That’s because the coordinates are now expressed in meters, representing positions along the X and Y axes of a 2D projected map. With this transformation complete, we’re ready to perform geospatial analysis.

Let’s start by calculating the distance between two sampling sites:

st_distance(sf_ft_utm)
## Units: [m]
##          1        2
## 1      0.0 106933.3
## 2 106933.3      0.0

These two sites are \(\sim\) 107 km away, shown in a matrix format.

2.5 How to choose CRS?

In a nutshell, geodetic CRS define how locations are mapped onto a 3D model of the Earth—typically an ellipsoid—preserving their original latitude and longitude values. These systems are ideal for storing and referencing geographic data globally but are not well suited for precise spatial measurements like area or distance due to the Earth’s curvature.

Projected CRS, on the other hand, transform geographic coordinates from the curved Earth onto a flat, 2D surface. Because no projection can preserve all spatial properties at once, the choice of a projected CRS should depend on the geometric attribute most important to your analysis—such as area, distance, or direction. Projected systems are also typically designed for use within specific regions, and applying them outside their intended coverage can introduce significant distortion. As a general guideline:

  • Use a UTM projection when working within a narrow longitudinal zone (e.g., for local studies).
  • Use a State Plane CRS for high-precision mapping within a single U.S. state.
  • Use a global equal-area projection (e.g., Mollweide or Albers Equal-Area)2 for continental or global-scale analyses.

Today, tools like generative AI can also serve as helpful guides in selecting an appropriate projected CRS by quickly narrowing down options based on your spatial extent, region, and analytical goals. Choosing the right projection helps ensure your spatial measurements and analyses are both accurate and interpretable.

2.6 Exercise

  1. Load the data (ref: Section 2.2)
    • Read the CSV file data_fish_deq.csv from the data folder using the readr::read_csv() function. This is a data set of fish fauna in North Carolina collected by the NC Department of Environmental Quality.
    • Assign the resulting tibble to an object called df_fish_deq.
    • Inspect df_fish_deq by printing and identify column names that contain unique site identifier and longitude/latitude information.
    • Create an object containing a column of unique site identifier along with longitude/latitude columns, and assign the resulting object to df_site_deq.
  2. Convert to an sf object (ref: Section 2.2)
    • Convert df_site_deq into an sf object using the sf::st_as_sf() function. Make sure to set the CRS to WGS84 (EPSG:4326) using the crs argument.
    • Assign the resulting spatial object to sf_site_deq.
    • Visualize sf_site_deq with mapview::mapview() function.
  3. Calculate the distance between two points (ref: Section 2.4)
    • Select the first two sites in sf_site_deq and assign it to sf_ft_deq.
    • Transform the CRS of sf_ft_deq to WGS 84 / UTM 17N using sf::st_transfrom(). Assign the resulting object to sf_ft_deq_utm.
    • Calculate the geographic distance between the two sites.
  4. Export the spatial object
    • Save the sf_site_deq object as an RDS file named sf_site_deq.rds in the data folder. Use saveRDS(sf_site_deq, file = "data/sf_site_deq.rds").