Map layers do not align with each other in R

Question

Recently, I have been asked to produce a map showing fault lines in Turkey after the unfortunate two earthquakes that took place today. Although I am not an expert in geospatial analysis, I have tried to produce this map with the following code chunk

library(elevatr)
library(terra)
library(tidyverse)
library(sf)
library(giscoR)
library(marmap)

# Getting the country's geospatial data

crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"

get_sf <- function(country_sf, country_transformed) {
  
  country_sf <- giscoR::gisco_get_countries(
    year = "2016",
    epsg = "4326",
    resolution = "10",
    country = "Turkey")
  
  country_transformed <- st_transform(country_sf, crs = crsLONGLAT)
  
  return(country_transformed)
}

country_transformed <- get_sf()

# Getting the fault lines data

shapename <- read_sf('gem_active_faults_harmonized.shp') %>%
  st_set_crs("+proj=longlat +datum=WGS84 +no_defs")

# Producing the map

ggplot() +
  geom_sf(data = country_transformed) +
  geom_sf(data = shapename, aes(color = "orange", size = 3)) +
  coord_sf(xlim = c(26, 45), ylim = c(36, 42)) +
  theme(legend.position = "none")

As can be seen from the code chunk, my workflow is based on 3 levels. First, getting the county's coordinates, and second getting the active fault lines data, which is available in this link and lastly merging these two datasets to produce the desired output. However, in the output produced as a result of this code, the layers do not seem align with each other. Therefore, do you have any suggestions to fix that problem? Thank you for your attention beforehand.

What's the CRS of your `gem_active_faults_harmonized.shp` file? You probably don't need to use `st_set_crs()` on it, because the shapefile should already contain the correct CRS info. If you need to change it to match the CRS of your other spatial layer, use `st_transform()`, not `st_set_crs()` — Stewart Macdonald, Feb 06 '23 at 22:53

Robert Hijmans · Accepted Answer · 2023-02-08T03:11:11.010

The map looks fine to me. Here is another way to make it.

Get the data

f <- 
"/vsicurl/https://github.com/GEMScienceTools/gem-global-active-faults/raw/master/geopackage/gem_active_faults_harmonized.gpkg"

library(terra)
flts <- vect(f)
wrld <- geodata::world(path=".")
wrld <- crop(wrld, ext(25, 45, 30, 43)+5)
epic <- cbind(37.203, 38.024)

With base-plot

plot(wrld, xlim=c(25, 45), ylim=c(30, 43), col=gray(0.95), back="azure", las=1)
lines(flts, col="red")
points(epic, pch=20, cex=2)

With ggplot

library(tidyterra)
library(ggplot2)
ggplot() +
  geom_spatvector(data = wrld) +
  geom_spatvector(data = flts, color = "red", size = 3) +
  geom_spatvector(data = vect(epic, crs=crs(flts))) +
  coord_sf(xlim = c(25, 45), ylim = c(30, 43)) +
  theme_bw() + theme(legend.position = "none")

And your approach with sf can be reduced to

library(sf)
country <- giscoR::gisco_get_countries(country = "Turkey")
faults <- sf::read_sf('gem_active_faults_harmonized.shp')
ggplot() +
  geom_sf(data = country) +
  geom_sf(data = faults, color = "orange", size = 3) +
  coord_sf(xlim = c(26, 45), ylim = c(36, 42)) +
  theme(legend.position = "none")

(using the shapefile, sf does not want to read the gpkg file).

Map layers do not align with each other in R

1 Answers1