Exploring the GeoThinned object

The output of thin_points() is a GeoThinned object, which contains: - x$retained: a list of logical vectors indicating which points were retained in each trial - x$method: the method used for thinning - x$params: the parameters used for thinning - x$original_data: the original dataset

You can apply basic methods for inspecting the object such as print(), summary(), and plot(). The summary() function provides detailed information about the thinning process, including the number of points retained in each trial, the nearest neighbor distance (NND), and the coverage of the retained points.

summary(quick_thin)
#> Summary of GeoThinneR Results
#> -----------------------------
#> Method used       : distance 
#> Distance metric   : haversine 
#> 
#> Number of points:
#>   Original              2000
#>   Thinned               1391
#>   Retention ratio      0.696
#> 
#> Nearest Neighbor Distances [km]
#>             Min 1st Qu. Median   Mean 3rd Qu.   Max
#> Original  0.650  11.011 17.159 18.216  24.265 61.83
#> Thinned  20.015  23.050 26.545 28.026  31.380 61.83
#> 
#> Spatial Coverage (Convex Hull Area) [km2]
#>   Original        2458025.112
#>   Thinned         2457093.024

We see that from 2000, we retained the 69.6% of points (1391), and the minimum distance between points has increased without reducing too much the spatial coverage or extent of the points compared to the original dataset.

plot(quick_thin)

If we check the number of points retained in each trial, we can see that the 5 trials return different numbers of points retained. This is due to the randomness of the algorithm when selecting which points to retain. You can use seed for reproducibility. By setting, all_trials = FALSE it will just return the trial that kept the most points.

# Number of kept points in each trial
sapply(quick_thin$retained, sum)
#> [1] 1387 1385 1388 1391 1387

You can access the thinned dataset that retained the maximum number of points using the largest() function.

head(largest(quick_thin))
#>          lon        lat  sp
#> 1 -4.2484496 -3.4032602 sp2
#> 2  5.7661027 -3.5548415 sp2
#> 3 -1.8204616 -3.5081961 sp2
#> 6 -9.0888700  1.1634277 sp1
#> 7  0.5621098 -0.5257711 sp1
#> 8  7.8483809 -4.4432328 sp1

If you want to access points from a specific trial, you can use the get_trial() function, which returns the thinned dataset for that trial.

head(get_trial(quick_thin, trial = 2)) 
#>          lon        lat  sp
#> 1 -4.2484496 -3.4032602 sp2
#> 2  5.7661027 -3.5548415 sp2
#> 3 -1.8204616 -3.5081961 sp2
#> 6 -9.0888700  1.1634277 sp1
#> 7  0.5621098 -0.5257711 sp1
#> 8  7.8483809 -4.4432328 sp1

Just in case you want to get the points that were removed from the original dataset, you can use the retained list to index the original data. The retained list contains logical vectors indicating which points were retained in each trial. You can use this to filter the original data and get the removed points.

trial <- 2
head(quick_thin$original_data[!quick_thin$retained[[trial]], ])
#>         lon         lat  sp
#> 4  7.660348  0.14434260 sp1
#> 5  8.809346 -0.07172694 sp1
#> 11 9.136667  3.50963224 sp2
#> 20 9.090073 -4.96464505 sp2
#> 21 7.790786  4.97424144 sp1
#> 32 8.045981  0.77635149 sp1

Also, you can use as_sf() to convert the thinned dataset to a sf object, which is useful for plotting and spatial analysis.

sf_thinned <- as_sf(quick_thin)
class(sf_thinned)
#> [1] "sf"         "data.frame"

Real-world example

GeoThinneR includes a dataset of loggerhead sea turtle (Caretta caretta) occurrence records in the Mediterranean Sea, sourced from GBIF and filtered. We now apply a 30 km thinning distance, returning only the largest trial (the one with the maximum number of points retained, all_trials set to FALSE).

# Apply spatial thinning to the real data
caretta_thin <- thin_points(
  data = caretta, # We will not specify lon_col, lat_col as they are in position 1 and 2
  lat_col = "decimalLatitude",
  lon_col = "decimalLongitude",
  method = "distance",
  thin_dist = 30, # Thinning distance in km,
  trials = 5,
  all_trials = FALSE,
  seed = 123
)

# Thinned dataframe stored in the first element of the retained list
identical(largest(caretta_thin), get_trial(caretta_thin, 1))
#> [1] TRUE

As you can see, with GeoThinneR, it’s very easy to apply spatial thinning to your dataset. In the next section we explain the different thinning methods and how to use them, and also the implemented optimizations for large-scale datasets.

Distance-based thinning

Distance-based thinning (method="distance") ensures that all retained points are separated by at least a user-defined minimum distance thin_dist (in km). For geographic coordinates, using distance = "haversine" (default) computes great-circle distances, or alternatively, the Euclidean distance between 3D Cartesian coordinates (depending on the method) to identify neighboring points within the distance threshold. The algorithm iterates through the dataset and removes points that are too close to each other.

In general, the choice of search method depends on the size of the dataset and the desired thinning distance:

• search_type="brute": Useful for small datasets, low scalability.
• search_type="kd_tree": More memory-efficient, better scalability for medium-sized datasets.
• search_type="local_kd_tree" (default): Faster than "kd_tree" and reduces memory usage. Can be parallelized to improve speed.
• search_type="k_estimation": Very fast when the number of points to remove is low or not very clustered, if not similar or worse than "local_kd_tree".

You can find the benchmarking results for each method in the package manuscript (see references).

dist_thin <- thin_points(sim_data, method = "distance", thin_dist = 25, trials = 3)
plot(dist_thin)

brute_thin <- thin_points(
  data = sim_data,
  method = "distance",
  search_type = "brute",
  thin_dist = 20,
  trials = 10,
  all_trials = FALSE,
  seed = 123
)
nrow(largest(brute_thin))
#> [1] 1391

kd_tree_thin <- thin_points(
  data = sim_data,
  method = "distance",
  search_type = "kd_tree",
  thin_dist = 20,
  trials = 10,
  all_trials = FALSE,
  seed = 123
)
nrow(largest(kd_tree_thin))
#> [1] 1391

local_kd_tree_thin <- thin_points(
  data = sim_data,
  method = "distance",
  search_type = "local_kd_tree",
  thin_dist = 20,
  trials = 10,
  all_trials = FALSE,
  seed = 123
)
nrow(largest(local_kd_tree_thin))
#> [1] 1391

k_estimation_thin <- thin_points(
  data = sim_data,
  method = "distance",
  search_type = "k_estimation",
  thin_dist = 20,
  trials = 10,
  all_trials = FALSE,
  seed = 123
)
nrow(largest(k_estimation_thin))
#> [1] 1391

Grid-based thinning

Grid sampling is a standard method where the area is divided into a grid, and n points are sampled from each grid cell. This method is very fast and memory-efficient. The only drawback is that you cannot ensure a minimum distance between points. There are two main ways to apply grid sampling in GeoThinneR: (i) define the characteristics of the grid, or (ii) pass your own grid as a raster (terra::SpatRaster).

For the first method, you can use the thin_dist parameter to define the grid cell size (the distance in km will be approximated to degrees to define the grid cell size, one degree roughly 111.32 km), or you can pass the resolution of the grid (e.g., resolution = 0.25 for 0.25x0.25-degree cells). If you want to align the grid with external data or covariate layers, you can pass the origin argument as a tuple of two values (e.g., c(0, 0)). Similarly, you can specify the coordinate reference system (CRS) of your grid (crs).

system.time(
grid_thin <- thin_points(
  data = sim_data,
  method = "grid",
  resolution = 1,
  origin = c(0, 0),
  crs = "epsg:4326",
  n = 1, # one point per grid cell
  trials = 10,
  all_trials = FALSE,
  seed = 123
))
#>    user  system elapsed 
#>    0.03    0.00    0.06
nrow(largest(grid_thin))
#> [1] 200

Alternatively, you can pass a terra::SpatRaster object, and that grid will be used for the thinning process.

rast_obj <- terra::rast(xmin = -10, xmax = 10, ymin = -5, ymax = 5, res = 1)
system.time(
grid_raster_thin <- thin_points(
  data = sim_data,
  method = "grid", 
  raster_obj = rast_obj,
  trials = 10,
  n = 1,
  all_trials = FALSE,
  seed = 123
))
#>    user  system elapsed 
#>    0.02    0.00    0.00
nrow(largest(grid_raster_thin))
#> [1] 200

Precision thinning

In this approach, coordinates are rounded to a certain precision to remove points that fall too close together. After removing points based on coordinate precision, the coordinate values are restored to their original locations. This is the simplest method and is more of a filtering when working with data from different sources with varying coordinate precision. To use it, you need to define the precision parameter, indicating the number of decimals to which the coordinates should be rounded.

system.time(
precision_thin <- thin_points(
  data = sim_data,
  method = "precision", 
  precision = 0,
  trials = 50,
  all_trials = FALSE,
  seed = 123
))
#>    user  system elapsed 
#>       0       0       0
nrow(largest(precision_thin))
#> [1] 230

These are the methods implemented in GeoThinneR. Depending on your specific dataset and research needs, one method may be more suitable than others.

Spatial thinning by group

In some cases, your dataset may include different groups, such as species, time periods, areas, or conditions, that you want to thin independently. The group_col parameter allows you to specify the column containing the grouping factor, and the thinning will be performed separately for each group. For example, in the simulated data where we have two species, we can use this parameter to thin each species independently:

all_thin <- thin_points(
  data = sim_data,
  thin_dist = 100,
  seed = 123
)
grouped_thin <- thin_points(
  data = sim_data,
  group_col = "sp",
  thin_dist = 100,
  seed = 123
)

nrow(largest(all_thin))
#> [1] 167
nrow(largest(grouped_thin))
#> [1] 321

Fixed number of points

What if you need to retain a fixed number of points that best covers the area where your data points are located, or for class balancing? The target_points parameter allows you to specify the number of points to keep, and the function will return that number of points spaced as separated as possible. Additionally, you can also set a thin_dist parameter so that no points closer than this distance will be retained. Currently, this approach is only implemented using the brute force method (method = "distance" and search_type = "brute"), so be cautious when applying it to very large datasets.

targeted_thin <- thin_points(
  data = caretta,
  lon_col = "decimalLongitude",
  lat_col = "decimalLatitude",
  target_points = 150,
  thin_dist = 30,
  all_trials = FALSE,
  seed = 123,
  verbose = TRUE
)
#> Starting spatial thinning at 2025-04-23 22:41:16 
#> Thinning using method: distance
#> For specific target points, brute force method is used.
#> Thinning process completed.
#> Total execution time: 3.57 seconds
nrow(largest(targeted_thin))
#> [1] 150

Select points by priority

In some scenarios, you may want to prioritize certain points based on a specific criterion, such as uncertainty, data quality, or any other variable. The priority parameter allows you to pass a vector representing the priority of each point. Currently, this feature can be used with the "grid" and "precision" methods.

For example, in the sea turtle data downloaded from GBIF, there is a column named coordinateUncertaintyInMeters. We can use this to prioritize points with lower uncertainty within each grid cell (in the grid method) or when rounding coordinates (in the precision method). Keep in mind that bigger uncertainty values represent less priority so we have to reverse this values.

grid_thin <- thin_points(
  data = caretta,
  lon_col = "decimalLongitude",
  lat_col = "decimalLatitude",
  method = "grid",
  resolution = 2,
  seed = 123
)

# Substracting the maximum - the highest uncertainty becomes the lowest priority and vice versa.
priority <- max(caretta$coordinateUncertaintyInMeters, na.rm = TRUE) - caretta$coordinateUncertaintyInMeters
priority_thin <- thin_points(
  data = caretta,
  lon_col = "decimalLongitude",
  lat_col = "decimalLatitude",
  method = "grid",
  resolution = 2,
  priority = priority,
  seed = 123
)

mean(largest(grid_thin)$coordinateUncertaintyInMeters, na.rm = TRUE)
#> [1] 35513.54
mean(largest(priority_thin )$coordinateUncertaintyInMeters, na.rm = TRUE)
#> [1] 17260.05