Ensuring Geospatial Accuracy: CRS Assignments and Transformations

Spatial Tutorial

This post aims to provide a comprehensive guide on correctly assigning coordinate reference systems (CRS) and handling spatial data with different CRS configurations. It will demonstrate how to project/transform geodata and their coordinates from one CRS to another, ensuring that all spatial objects (maps, rasters, shapefiles) align accurately in a common reference frame.

Moritz Wenzler-Meya (IZW Berlin)https://ecodynizw.github.io/EcoDynIZW
2025-02-28

Load packages

Before starting, we need to load the necessary R packages. It is essential to load only the required packages to avoid conflicts, especially between packages like terra and raster, which share function names.

library(sf)
library(terra)

Set paths

Next, we define the paths to our working directories. Adjust these paths according to your folder structure.

#main_dir <- ""

Set the paths for the necessary datasets from our d6_teaching_collection github repository (https://github.com/EcoDynIZW/d6_teaching_collection). When naming variables, it is best practice to avoid capital letters and use short, descriptive names for ease of typing and readability.

subdir_berlin <- paste(main_dir, "data_berlin", sep = "/")
subdir_tif <- paste(subdir_berlin, "geo_raster_current_gtif", sep = "/")
subdir_asc <- paste(subdir_berlin, "geo_raster_current_asc", sep = "/")
subdir_animals <- paste(subdir_berlin, "animal_data", sep = "/")
subdir_shape <- paste(main_dir, "data_borneo", "geo_vector", sep = "/")

Load data

We will load different data types and demonstrate how to assign and transform the coordinate reference systems.

Load csv files

# excel/csv files
boar <- read.csv(paste(subdir_animals, 
                       "data_wb_melden_en.csv", 
                       sep = "/"))
head(boar)
  observation id    gps_y    gps_x                 date year
1           1  9 52.64371 13.31417 2013-06-07T00:00:00Z 2013
2           2 11 52.61391 13.22447 2013-06-08T00:00:00Z 2013
3           3 26 52.44837 13.69075 2013-07-24T00:00:00Z 2013
4           4 32 52.46341 13.24752 2013-08-08T00:00:00Z 2013
5           5 56 52.61388 13.22801 2013-08-15T00:00:00Z 2013
6           6 85 52.38521 13.63526 2013-10-07T00:00:00Z 2013
  group_size n_piglets weather
1          5         2  cloudy
2          5        10   sunny
3         10         4   sunny
4          1         1   sunny
5         10         4   sunny
6         10         7  cloudy

Create vector data from a dataframe

If our spatial data is stored in a dataframe, we can convert it into an sf object using st_as_sf, specifying the coordinate columns and the CRS.

boar_sf <- st_as_sf(x = boar,
                    coords = c("gps_x", "gps_y"),
                    crs = 4326, # epsg code to use for crs
                    sf_column_name = "geometry", # standard name
                    remove = FALSE) # if you want to keep the coordinate columns

Load vector data

Now, we load vector data. The st_read function from the sf package allows us to read and inspect these files.

borneo_sf <- st_read(paste(subdir_shape ,"borneo_admin.shp", sep = "/"))
Reading layer `borneo_admin' from data source 
  `C:\Users\wenzler\Documents\GitHub\d6_teaching_collection\data\data_borneo\geo_vector\borneo_admin.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 158 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 108.5986 ymin: -4.943054 xmax: 119.2697 ymax: 7.380556
Geodetic CRS:  WGS 84
borneo_sf
Simple feature collection with 158 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 108.5986 ymin: -4.943054 xmax: 119.2697 ymax: 7.380556
Geodetic CRS:  WGS 84
First 10 features:
   ID_0 ISO   NAME_0 ID_1 NAME_1  ID_2     NAME_2 VARNAME_2 NL_NAME_2
1   158 MYS Malaysia 1988  Sabah 23051 Lahad Datu      <NA>      <NA>
2   158 MYS Malaysia 1988  Sabah 23052      Papar      <NA>      <NA>
3   158 MYS Malaysia 1988  Sabah 23053  Penampang      <NA>      <NA>
4   158 MYS Malaysia 1988  Sabah 23054 Pensiangan      <NA>      <NA>
5   158 MYS Malaysia 1988  Sabah 23055      Pitas      <NA>      <NA>
6   158 MYS Malaysia 1988  Sabah 23056      Ranau      <NA>      <NA>
7   158 MYS Malaysia 1988  Sabah 23057   Sandakan      <NA>      <NA>
8   158 MYS Malaysia 1988  Sabah 23058   Semporna      <NA>      <NA>
9   158 MYS Malaysia 1988  Sabah 23059   Sipitang      <NA>      <NA>
10  158 MYS Malaysia 1988  Sabah 23060   Tambunan      <NA>      <NA>
     HASC_2 CC_2   TYPE_2 ENGTYPE_2 VALIDFR_2 VALIDTO_2 REMARKS_2
1  MY.SA.LD <NA> District  District   Unknown   Present      <NA>
2  MY.SA.PA <NA> District  District   Unknown   Present      <NA>
3  MY.SA.PE <NA> District  District   Unknown   Present      <NA>
4      <NA> <NA> District  District   Unknown   Unknown      <NA>
5  MY.SA.PI <NA> District  District   Unknown   Present      <NA>
6  MY.SA.RA <NA> District  District   Unknown   Present      <NA>
7  MY.SA.SA <NA> District  District   Unknown   Present      <NA>
8  MY.SA.SE <NA> District  District   Unknown   Present      <NA>
9  MY.SA.SI <NA> District  District   Unknown   Present      <NA>
10 MY.SA.TM <NA> District  District   Unknown   Present      <NA>
   Shape_Leng Shape_Area                       geometry
1   6.7475367 0.62106012 MULTIPOLYGON (((118.2299 4....
2   1.6727672 0.10065591 MULTIPOLYGON (((116.0144 5....
3   0.9548857 0.03951545 MULTIPOLYGON (((116.0477 5....
4   4.0257859 0.49725547 MULTIPOLYGON (((116.5357 5....
5   2.9223456 0.11955352 MULTIPOLYGON (((117.2762 6....
6   2.4815704 0.24027591 MULTIPOLYGON (((116.9182 6....
7   7.1555231 0.18354676 MULTIPOLYGON (((117.9468 5....
8   6.0198462 0.09494946 MULTIPOLYGON (((118.6294 4....
9   2.9190984 0.22604467 MULTIPOLYGON (((115.7132 5....
10  1.5307563 0.12152936 MULTIPOLYGON (((116.6079 5....

Load GeoTIFF files

GeoTIFF files conveniently store CRS information within the file itself. We can extract the CRS specification from one of our GeoTIFFs.

temp <- rast(paste0(subdir_tif, "/summer_temp_100m_3035.tif"))
temp
class       : SpatRaster 
dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source      : summer_temp_100m_3035.tif 
name        : summer_temp_100m_3035 
min value   :                  2052 
max value   :                  3133

Load ASC files

Here, we load an ASC raster file. Unlike GeoTIFFs, ASC files often lack embedded CRS information, so we need to define it manually.

tree_cover <- rast(paste0(subdir_asc, "/tree_cover_density_2012_100m_3035.asc"))
tree_cover
class       : SpatRaster 
dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
coord. ref. :  
source      : tree_cover_density_2012_100m_3035.asc 
name        : tree_cover_density_2012_100m_3035

What is missing?

Set the coordinate reference system for an ASC file

If a CRS is missing, we must define it based on prior knowledge of the dataset. By checking the coordinate extent, we can infer that these are not latitude and longitude coordinates but belong to the EPSG:3035 projection.

You can check here what it stands for: https://epsg.io/3035

This website also provides the full-text specification (scroll down to ‘export’ and then klick on ‘PROJ.4’):

+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs But his method is quite old and needs a lot of characters what can leads to more errors.

crs(tree_cover) <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80                                       +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

A better option is to use a crs from a file you already using.

# or 
crs(tree_cover) <- crs(temp)

or using the state of the art wkt-crs (Well-known text representation of coordinate reference systems).

# or
crs(tree_cover) <- st_crs(3035)$wkt

or staying in the terra envirinment using just the epsg code.

# or 
crs(tree_cover) <- crs("epsg:3035")

After assigning the CRS, we verify if it has been correctly applied.

tree_cover
class       : SpatRaster 
dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
source      : tree_cover_density_2012_100m_3035.asc 
name        : tree_cover_density_2012_100m_3035

How to project / transform data

When working with spatial data from different sources, their coordinate systems must match to ensure accurate spatial overlays. Below, we explore how to transform different types of spatial data to a common CRS.

Project an sf object (vector data)

Transforming vector data is more reliable than raster transformations. Below, we reproject an sf object to a new CRS.

borneo_sf_3035 <- st_transform(borneo_sf, 3035)
borneo_sf_3035
Simple feature collection with 158 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 13226770 ymin: 4169872 xmax: 14197800 ymax: 6257863
Projected CRS: ETRS89-extended / LAEA Europe
First 10 features:
   ID_0 ISO   NAME_0 ID_1 NAME_1  ID_2     NAME_2 VARNAME_2 NL_NAME_2
1   158 MYS Malaysia 1988  Sabah 23051 Lahad Datu      <NA>      <NA>
2   158 MYS Malaysia 1988  Sabah 23052      Papar      <NA>      <NA>
3   158 MYS Malaysia 1988  Sabah 23053  Penampang      <NA>      <NA>
4   158 MYS Malaysia 1988  Sabah 23054 Pensiangan      <NA>      <NA>
5   158 MYS Malaysia 1988  Sabah 23055      Pitas      <NA>      <NA>
6   158 MYS Malaysia 1988  Sabah 23056      Ranau      <NA>      <NA>
7   158 MYS Malaysia 1988  Sabah 23057   Sandakan      <NA>      <NA>
8   158 MYS Malaysia 1988  Sabah 23058   Semporna      <NA>      <NA>
9   158 MYS Malaysia 1988  Sabah 23059   Sipitang      <NA>      <NA>
10  158 MYS Malaysia 1988  Sabah 23060   Tambunan      <NA>      <NA>
     HASC_2 CC_2   TYPE_2 ENGTYPE_2 VALIDFR_2 VALIDTO_2 REMARKS_2
1  MY.SA.LD <NA> District  District   Unknown   Present      <NA>
2  MY.SA.PA <NA> District  District   Unknown   Present      <NA>
3  MY.SA.PE <NA> District  District   Unknown   Present      <NA>
4      <NA> <NA> District  District   Unknown   Unknown      <NA>
5  MY.SA.PI <NA> District  District   Unknown   Present      <NA>
6  MY.SA.RA <NA> District  District   Unknown   Present      <NA>
7  MY.SA.SA <NA> District  District   Unknown   Present      <NA>
8  MY.SA.SE <NA> District  District   Unknown   Present      <NA>
9  MY.SA.SI <NA> District  District   Unknown   Present      <NA>
10 MY.SA.TM <NA> District  District   Unknown   Present      <NA>
   Shape_Leng Shape_Area                       geometry
1   6.7475367 0.62106012 MULTIPOLYGON (((13458143 60...
2   1.6727672 0.10065591 MULTIPOLYGON (((13355721 58...
3   0.9548857 0.03951545 MULTIPOLYGON (((13352723 58...
4   4.0257859 0.49725547 MULTIPOLYGON (((13424936 58...
5   2.9223456 0.11955352 MULTIPOLYGON (((13283399 60...
6   2.4815704 0.24027591 MULTIPOLYGON (((13325784 59...
7   7.1555231 0.18354676 MULTIPOLYGON (((13376142 60...
8   6.0198462 0.09494946 MULTIPOLYGON (((13511000 60...
9   2.9190984 0.22604467 MULTIPOLYGON (((13418281 57...
10  1.5307563 0.12152936 MULTIPOLYGON (((13359817 59...
#synonyms:
#borneo_sf_3035 <- st_transform(borneo_sf, crs(tree_cover))
#borneo_sf_3035 <- st_transform(borneo_sf, st_crs(3035)$wkt)

Project a raster

It is generally advisable to transform vector data instead of raster data to preserve spatial accuracy. More specific, it can lead to unusual cell sizes which can can corrupt your calculations with these files.

However, for demonstration purposes, we apply a CRS transformation to a raster.

In this case we want to change the CRS from the temperature raster from EPSG 3035 to EPSG 25833 (UTM).

temp_utm <- project(x = temp, y = st_crs(25833)$wkt) 
temp_utm
class       : SpatRaster 
dimensions  : 614, 695, 1  (nrow, ncol, nlyr)
resolution  : 99.96378, 99.96378  (x, y)
extent      : 358826.2, 428301, 5788287, 5849665  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 33N (EPSG:25833) 
source(s)   : memory
name        : summer_temp_100m_3035 
min value   :              2066.307 
max value   :              3133.000

Handling raster transformations

If transforming a raster is unavoidable, using the resample function afterward can help maintain cell alignment and structure. Nevertheless, you need a raster with the correct cellsize and projection as a template for the resampling prozess.

temp_3035 <- project(x = temp_utm, y = st_crs(3035)$wkt) 
temp_3035_res <- resample(x = temp_3035, y = temp, method = "near") # change method depending on the data input
temp_3035_res
class       : SpatRaster 
dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
source(s)   : memory
varname     : summer_temp_100m_3035 
name        : summer_temp_100m_3035 
min value   :              2073.621 
max value   :              3133.000 
temp
class       : SpatRaster 
dimensions  : 570, 657, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 4521040, 4586740, 3243800, 3300800  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source      : summer_temp_100m_3035.tif 
name        : summer_temp_100m_3035 
min value   :                  2052 
max value   :                  3133

After resampling, some values may have changed due to interpolation. Always check the transformed raster carefully to ensure its usability.

Ensuring that all spatial data align within a common coordinate reference system is crucial for accurate spatial analysis. By correctly assigning, transforming, and verifying CRS for different data formats, we can prevent misalignment issues and maintain spatial integrity in geospatial projects.

Citation

For attribution, please cite this work as

Wenzler-Meya (2025, Feb. 28). Ecological Dynamics: Ensuring Geospatial Accuracy: CRS Assignments and Transformations. Retrieved from https://ecodynizw.github.io/posts/common-mistakes-with-geodata/

BibTeX citation

@misc{wenzler-meya2025ensuring,
  author = {Wenzler-Meya, Moritz},
  title = {Ecological Dynamics: Ensuring Geospatial Accuracy: CRS Assignments and Transformations},
  url = {https://ecodynizw.github.io/posts/common-mistakes-with-geodata/},
  year = {2025}
}

© Department of Ecological Dynamics, Leibniz Institute for Zoo- and Wildlife Research 2020—2025. All rights reserved.

Impressum