data("quakes")
+crs <- st_crs("EPSG:3460")
-quakes_sf <- st_as_sf(quakes, coords = c("long", "lat"), crs = "WGS84")
+quakes_sf <- st_as_sf(quakes, coords = c("long", "lat"), crs = "WGS84") |>
+ st_transform(crs = crs)
quake_map <- ggplot() +
geom_sf(aes(color = depth, size = mag), data = quakes_sf) +
- scale_color_viridis_c() +
scale_size_binned() +
- theme_classic() +
- coord_sf(xlim = c(160, 195))
+ theme_classic()
quake_map
+> Transform your dataset into a projected cooridnate reference system
+to deal with units such as m, km, miles etc. Otherwise the units will be
+in degrees.
profile_pts <- data.frame(lon = c(160, -170), lat = c(-15, -24)) |>
st_as_sf(coords = c("lon", "lat"), crs = "WGS84") |>
- st_shift_longitude()Combine the two points to a line and add the profile line to the map:
profile_l <- profile_line(profile_pts)
quake_map +
- geom_sf(data = profile_l, lty = 2)
-#> Coordinate system already present. Adding new coordinate system, which will
-#> replace the existing one.
profile_azimuth(profile_l)
#> 112.3832 [°]
profile_length(profile_l)
-#> 3291737 [m]
data.frame(lon = 160, lat = 15) |>
st_as_sf(coords = c("lon", "lat"), crs = "WGS84") |>
- st_shift_longitude() |>
- profile_points(profile.azimuth = 112, profile.length = set_units(50, degree))
+ st_transform(crs = crs) |>
+ profile_points(profile.azimuth = 112, profile.length = set_units(8000, km))
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: POINT
#> Dimension: XY
-#> Bounding box: xmin: 160 ymin: -31.35919 xmax: 178.7303 ymax: 15
-#> Geodetic CRS: WGS 84
-#> geometry
-#> 1 POINT (160 15)
-#> 2 POINT (178.7303 -31.35919)Note that the unit of
profile.lengthdepends on the coordinate reference system and must be in degree for lon-lat @@ -161,9 +162,12 @@Determine the dist plotting these axes against each other::
+ geom_point() + + scale_x_continuous(breaks = seq(0, 3000, 250)) + + scale_y_continuous(breaks = seq(-3000, 3000, 250)) + + coord_fixed()quakes_profile |> - ggplot(aes(X, Y, color = depth)) + + ggplot(aes(X / 1000, Y / 1000, color = depth)) + geom_hline(yintercept = 0) + - geom_point()@@ -172,8 +176,8 @@
In this plot the profile line is a horizontal line (
X=0).Determine the dist
quakes_profile_filtered <- filter( quakes_profile, - abs(Y) <= 10, - X <= 7 + abs(Y) <= 750 * 1000, + X >= 1500 * 1000 )@@ -186,7 +190,7 @@+ labs(x = "Distance along profile (m)", y = "Depth (km)", color = "Depth (km)")Plot data along profile scale_size_binned("Richter Magnitude") + scale_y_reverse() + scale_x_continuous(guide = guide_axis(position = "top")) + - labs(x = "Distance along profile (degree)", y = "Depth (km)", color = "Depth (km)")
One way to show the distance from the profile in the same plot is by controlling the size (and/or the opacity) of the points. Here I am using @@ -200,7 +204,7 @@
Plot data along profile scale_alpha_continuous("Distance from profile (degree)", range = c(1, .1)) + scale_y_reverse() + scale_x_continuous(guide = guide_axis(position = "top")) + - labs(x = "Distance along profile (degree)", y = "Depth (km)") + labs(x = "Distance along profile (m)", y = "Depth (km)")
