inst/doc/geomattR.R

In geomattR: Calculate Geometric Attributes of Spatial Polygons

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE
)

## ----setup--------------------------------------------------------------------
# # Install if you haven't already
# # install.packages("devtools")
# # devtools::install_github("gortegasolis/geomattR")
# 
# library(geomattR)
# library(terra)

## ----quick_start--------------------------------------------------------------
# # Create a simple polygon (1x1 degree square at equator)
# coords <- cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))
# polygon <- vect(coords, type = "polygon", crs = "EPSG:4326")
# 
# # Calculate all available metrics
# result <- calculate_geometric_attributes(polygon)
# 
# # View the results
# print(result)
# 
# # Or view specific columns
# print(result[, c("area", "perimeter", "compactness")])

## ----size_metrics-------------------------------------------------------------
# # Calculate area and perimeter
# polygon_with_size <- calculate_geometric_attributes(
#   polygon,
#   metrics = c("area", "perimeter", "hole_area", "hole_area_pct")
# )
# 
# # Area in square meters (using geodesic calculations)
# print(polygon_with_size$area)
# 
# # Perimeter in meters
# print(polygon_with_size$perimeter)

## ----shape_metrics------------------------------------------------------------
# # Calculate shape characteristics
# polygon_with_shape <- calculate_geometric_attributes(
#   polygon,
#   metrics = c(
#     "compactness", # Polsby-Popper (0-1, 1=perfect circle)
#     "reock", # Reock compactness
#     "shape_index", # Dimensionless shape complexity
#     "circularity_ratio", # Based on max length
#     "fractaldimension", # Boundary complexity
#     "elongation_rectangle" # mean(long sides) / mean(short sides) from the minimum bounding rectangle
#   )
# )
# 
# # Perfect circle has compactness = 1
# # Irregular shapes have lower values
# print(polygon_with_shape[, c("compactness", "shape_index")])

## ----orientation_metrics------------------------------------------------------
# # Calculate orientation
# polygon_with_orient <- calculate_geometric_attributes(
#   polygon,
#   metrics = c("bearing", "northerness", "ew_length", "ns_length", "maxlength")
# )
# 
# # Bearing: axial orientation of the maximum length line, always measured
# # south to north, so it lies within [-90, 90] (0 = N-S axis, -90/90 = E-W axis)
# print(polygon_with_orient$bearing)
# 
# # Northerness: cosine of bearing (0 = E-W axis, 1 = N-S axis)
# print(polygon_with_orient$northerness)
# 
# # Extents in meters: maximum distances from East to West and North to South
# print(paste("EW Extent:", round(polygon_with_orient$ew_length, 0), "m"))
# print(paste("NS Extent:", round(polygon_with_orient$ns_length, 0), "m"))

## ----complexity_metrics-------------------------------------------------------
# polygon_with_complexity <- calculate_geometric_attributes(
#   polygon,
#   metrics = c("sinuosity", "fractaldimension")
# )
# 
# # Sinuosity: perimeter / max_length (>1 indicates rougher boundary)
# print(polygon_with_complexity$sinuosity)
# 
# # Fractal dimension (FRAGSTATS convention): boundary complexity, approaching
# # 1 for simple square-like shapes and 2 for highly convoluted boundaries
# print(polygon_with_complexity$fractaldimension)

## ----geometry_metrics---------------------------------------------------------
# # Get centroid coordinates and structure info
# polygon_with_geom <- calculate_geometric_attributes(
#   polygon,
#   metrics = c(
#     "decimallongitude",
#     "decimallatitude",
#     "num_holes",
#     "num_polygons"
#   )
# )
# 
# print(polygon_with_geom[, c("decimallongitude", "decimallatitude")])

## ----multiple_polygons--------------------------------------------------------
# # Create multiple polygons
# coords1 <- cbind(c(0, 0, 1, 1, 0), c(0, 1, 1, 0, 0))
# coords2 <- cbind(c(1, 1, 2, 2, 1), c(1, 2, 2, 1, 1))
# 
# poly1 <- vect(coords1, type = "polygon", crs = "EPSG:4326")
# poly2 <- vect(coords2, type = "polygon", crs = "EPSG:4326")
# 
# # Combine into single SpatVector
# polygons <- rbind(poly1, poly2)
# 
# # Calculate metrics for all polygons at once
# results <- calculate_geometric_attributes(
#   polygons,
#   metrics = c("area", "perimeter", "compactness")
# )
# 
# # Results have same row count as input
# print(nrow(results)) # 2
# 
# # All metrics calculated
# print(results[, c("area", "perimeter", "compactness")])

## ----parallel, eval=FALSE-----------------------------------------------------
# library(parallel)
# 
# # Create multiple polygons
# polygons <- rbind(poly1, poly2) # From previous example
# 
# # Set up cluster (two workers here; adjust to your machine's resources)
# cl <- makeCluster(2)
# 
# # Process in parallel using the main function
# results <- calculate_geometric_attributes(
#   polygons,
#   metrics = "all",
#   cl = cl
# )
# 
# # Clean up
# stopCluster(cl)
# 
# # Results should be equivalent to sequential processing
# # (allowing for tiny floating-point differences)
# print(head(results))

## ----crs_handling-------------------------------------------------------------
# # Create polygon in geographic CRS (WGS84/EPSG:4326)
# poly_geographic <- vect(coords1, type = "polygon", crs = "EPSG:4326")
# 
# # Create same polygon in projected CRS (UTM Zone 30N)
# poly_projected <- project(poly_geographic, "EPSG:32630")
# 
# # Calculate metrics on both
# result_geo <- calculate_geometric_attributes(poly_geographic, metrics = "area")
# result_proj <- calculate_geometric_attributes(poly_projected, metrics = "area")
# 
# # Both should give similar results (using geodesic calculations by default)
# print(paste("Geographic CRS area:", round(result_geo$area, 0), "m²"))
# print(paste("Projected CRS area:", round(result_proj$area, 0), "m²"))
# 
# # Difference should be minimal due to default geodesic calculations
# print(paste(
#   "Difference:",
#   round(abs(result_geo$area - result_proj$area), 0),
#   "m²"
# ))