diff --git a/articles/geoprofiler.html b/articles/geoprofiler.html index f7dd98e..59c3ea1 100644 --- a/articles/geoprofiler.html +++ b/articles/geoprofiler.html @@ -69,7 +69,7 @@

2024-12-04

library(dplyr) theme_set(theme_classic()) -options(ggplot2.continuous.colour="viridis") +options(ggplot2.continuous.colour = "viridis") options(ggplot2.continuous.fill = "viridis")

Load example data @@ -79,18 +79,21 @@

Load example data
 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.

Define a profile @@ -105,22 +108,20 @@

Profile from two known points
 profile_pts <- data.frame(lon = c(160, -170), lat = c(-15, -24)) |>
   st_as_sf(coords = c("lon", "lat"), crs = "WGS84") |>
-  st_shift_longitude()

+ st_transform(crs = crs)

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.
+ geom_sf(data = profile_l, lty = 2) 

 profile_azimuth(profile_l)
 #> 112.3832 [°]
 profile_length(profile_l)
-#> 3291737 [m]
+#> 3336889 [m]

Profile from direction and length from one point @@ -130,16 +131,16 @@

Profile from direction 
 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)
+#> Bounding box: xmin: -48308.26 ymin: 210366.1 xmax: 2948544 ymax: 7627837 +#> Projected CRS: Fiji 1986 / Fiji Map Grid +#> geometry +#> 1 POINT (-48308.26 7627837) +#> 2 POINT (2948544 210366.1)

Note that the unit of profile.length depends 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::

 quakes_profile |>
-  ggplot(aes(X, Y, color = depth)) +
+  ggplot(aes(X / 1000, Y / 1000, color = depth)) +
   geom_hline(yintercept = 0) +
-  geom_point()
+ geom_point() + + scale_x_continuous(breaks = seq(0, 3000, 250)) + + scale_y_continuous(breaks = seq(-3000, 3000, 250)) + + coord_fixed()

In this plot the profile line is a horizontal line (X=0).

@@ -172,8 +176,8 @@

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 @@

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

+ labs(x = "Distance along profile (m)", 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)")