Search code examples
rgeospatialr-sfgeom-sf

Map layers do not align with each other in R

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.

enter image description here

Solution

  • 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)

    enter image description here

    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")

    enter image description here

    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).