[WIP] Proposed vector data format and application to LandIQ data

modules/data.land/R/landiq2std.R

@@ -0,0 +1,128 @@
+#' Convert a LandIQ Shapefile into a Standardized Format consisting of
+#' a GeoPackage file to store geospatial information and an associated
+#' CSV file with attributes
+#'
+#' This function reads a LandIQ crop map shapefile downloaded from
+#' https://www.landiq.com/land-use-mapping and processes the data into
+#' a standardized GeoPackage and CSV format.
+#'
+#' @param input_file Character. Path to the input Shapefile
+#' (GeoPackage created by shp2gpkg is also valid)
+#' @param output_gpkg Character. Path to the output GeoPackage
+#' @param output_csv Character. Path to the output CSV.
+#'
+#' @return Invisibly returns a list with paths to the output files.
+#' @details
+#' This function:
+#' - Reads a Shapefile using the `sf` package.
+#' - Calculates centroids to extract latitude and longitude.
+#' - Converts `Acres` column to hectares (`ha`).
+#' - Extracts crop information and assigns a PFT (Plant Functional Type).
+#' - Outputs a standardized GeoPackage with columns `id`, `geometry`, `lat`, `lon`, and `area_ha`.
+#' - Outputs a CSV with columns `year`, `crop`, `pft`, `source`, and `notes`.
+#'
+#' Note: TODO provide crop-->PFT mapping as an external file using either
+#' defaults stored with the PEcAn.data.land package or an optional custom
+#' mapping
+#'
+#' @examples
+#' # Specify input and output file paths
+#' input_file <- "path/to/landiq.shp"
+#' output_gpkg <- "path/to/your_output.gpkg"
+#' output_csv <- "path/to/your_output.csv"
+#'
+#' # Run the conversion function
+#' landiq2std(input_file, output_gpkg, output_csv)
+#'
+#' @importFrom sf st_read st_transform st_centroid st_coordinates st_write
+#' @importFrom dplyr mutate select rename filter case_when distinct
+#' @export
+landiq2std <- function(input_file, output_gpkg, output_csv) {
+  # Check input file format
+  if (grepl(pattern = "\\.shp$", input_file)) {
+    PEcAn.logger::logger.info(
+      "Converting Shapefile:", basename(input_file),
+      "to GeoPackage:", basename(output_gpkg)
+    )
+    # read in, repair geometries, write out repaired geopackage
+    shp2gpkg(input_file, output_gpkg, overwrite = TRUE)
+    input_file <- output_gpkg # now gpkg is input file
+  }
+
+  # Read the Shapefile
+  landiq_polygons <- sf::st_read(input_file, quiet = FALSE)
+
+  # Determine crop year from column name
+  crop_col <- grep("Crop[0-9]{4}", colnames(landiq_polygons), value = TRUE)
+  year <- gsub("Crop", "", crop_col)
+  # Check required columns
+  required_cols <- c("Acres", crop_col, "Source", "Comments", "County", "geom")
+  missing_cols <- setdiff(required_cols, colnames(landiq_polygons))
+  if (length(missing_cols) > 0) {
+    stop("Input file is missing the following columns: ", paste(missing_cols, collapse = ", "))
+  }
+
+  # possible spedup by pre-computing lat and lon
+  landiq_polygons_updated <- landiq_polygons |>
+    sf::st_transform(4326) |>
+    dplyr::mutate(
+      year = year,
+      lon = sf::st_coordinates(sf::st_centroid(geom))[, "X"],
+      lat = sf::st_coordinates(sf::st_centroid(geom))[, "Y"],
+      area_ha = PEcAn.utils::ud_convert(Acres, "acre", "ha")
+    ) |>
+    dplyr::rowwise() |>
+    dplyr::mutate( # generate ids rowwise separately because
+      # rowwise geospatial operations are very slow
+      id = digest::digest(geom, algo = "xxhash64")
+    ) |>
+    dplyr::rename(county = County)
+
+  # Process data for GeoPackage
+  gpkg_data <- landiq_polygons_updated |>
+    dplyr::select(id, geom, lat, lon, area_ha, county)
+
+  # Process data for CSV
+  csv_data <- landiq_polygons_updated |>
+    tidyr::as_tibble() |>
+    dplyr::mutate(
+      crop = .data[[crop_col]],
+      pft = case_when(
+        crop %in% c(
+          "Cherries", "Almonds", "Plums, Prunes and Apricots",
+          "Walnuts", "Citrus", "Miscellaneous Deciduous", "Pears", "Olives",
+          "Apples", "Pistachios", "Bush Berries", "Peaches/Nectarines",
+          "Miscellaneous Subtropical Fruits", "Pomegranates"
+        ) ~ "woody perennial crop",
+        TRUE ~ paste0("no PFT for ", crop)
+      )
+    ) |>
+    dplyr::rename(
+      source = Source,
+      notes = Comments
+    ) |>
+    dplyr::select(id, lat, lon, year, crop, pft, source, notes)
+
+  # Warn about crops without a PFT
+  unassigned_pft <- csv_data |>
+    filter(grepl("no PFT for", pft)) |>
+    distinct(crop, pft)
+  if (nrow(unassigned_pft) > 0) {
+    PEcAn.logger::logger.warn( 
+      "The following crops do not have a PFT assigned:",
+      paste(unassigned_pft$crop, collapse = ", ")
+    )
+  }
+
+  # Write outputs
+  file.remove(output_gpkg, output_csv)
+  sf::st_write(gpkg_data,
+    output_gpkg,
+    layer = "sites", # analogous to BETYdb table name
+    quiet = FALSE
+  )
+  readr::write_csv(csv_data, output_csv)
+
+  # Return paths to output files
+  invisible(list(GeoPackage = output_gpkg, CSV = output_csv))
+}

modules/data.land/R/shp2gpkg.R

@@ -0,0 +1,84 @@
+#' Convert Shapefile to GeoPackage
+#'
+#' This function converts a Shapefile to a GeoPackage using the `sf` package.
+#'
+#' @param input_shp Character. Path to the input Shapefile (e.g., `"data/myfile.shp"`).
+#' @param output_gpkg Character. Path to the output GeoPackage (e.g., `"output/myfile.gpkg"`).
+#' @param layer_name Character. Optional. Name of the layer in the GeoPackage. Defaults to the base name of the Shapefile.
+#' @param overwrite Logical. Whether to overwrite an existing layer in the GeoPackage. Defaults to `TRUE`.
+#'
+#' @return Invisibly returns the path to the GeoPackage (`output_gpkg`) upon successful conversion.
+#' @details
+#' This function reads a Shapefile, converts it to an `sf` object, and writes it to a GeoPackage.
+#' The `sf` package handles the conversion and ensures spatial data integrity.
+#'
+#' @examples
+#' # Convert 'roads.shp' to 'roads.gpkg' with the default layer name
+#' shp2gpkg("data/roads.shp", "output/roads.gpkg")
+#'
+#' # Convert with a custom layer name
+#' shp2gpkg("data/roads.shp", "output/roads.gpkg", layer_name = "custom_layer")
+#'
+#' # Prevent overwriting existing layers
+#' shp2gpkg("data/roads.shp", "output/roads.gpkg", overwrite = FALSE)
+#'
+#' @importFrom sf st_read st_write st_is_valid st_make_valid
+#' @export
+shp2gpkg <- function(input_shp, output_gpkg, layer_name = NULL, overwrite = TRUE) {
+  # Load the Shapefile
+  shapefile <- sf::st_read(input_shp, quiet = TRUE)
+
+  # Check validity of geometries
+  is_geometry_valid <- sf::st_is_valid(shapefile)
+  total_geometries <- length(is_geometry_valid)
+  n_invalid <- sum(!is_geometry_valid)
+  percent_invalid <- (n_invalid / total_geometries) * 100
+
+  # Log invalid geometry statistics
+
+  PEcAn.logger::logger.info("Total records in shapefile: ", total_geometries)
+
+  if (n_invalid > 0) {
+    # Report invalid geometry statistics
+    PEcAn.logger::logger.info(
+      n_invalid,
+      "(", round(percent_invalid, 2), "%)",
+      "of geometries are invalid"
+    )
+    PEcAn.logger::logger.info(". Attempting to repair invalid geometries")
+
+    # Attempt to repair invalid geometries
+    repaired_shapefile <- sf::st_make_valid(shapefile)
+    is_repair_successful <- sf::st_is_valid(repaired_shapefile)
+    n_repaired <- n_invalid - sum(!is_repair_successful)
+    n_remaining_invalid <- sum(!is_repair_successful)
+
+    # Report repair results
+    message("Repaired ", n_repaired, " geometries successfully.")
+    if (n_remaining_invalid > 0) {
+      PEcAn.logger::logger.warn("Could not repair ", n_remaining_invalid, " geometries. They will be excluded.")
+    } else {
+      PEcAn.logger::logger.info("All geometries were repaired successfully.")
+    }
+
+    # Retain only valid geometries after repair
+    shapefile <- repaired_shapefile[is_repair_successful, ]
+  }
+
+  # Determine layer name from input file if not provided
+  if (is.null(layer_name)) {
+    layer_name <- tools::file_path_sans_ext(basename(input_shp))
+  }
+
+  # Write to GeoPackage
+  sf::st_write(
+    shapefile,
+    output_gpkg,
+    layer = layer_name,
+    delete_layer = overwrite,
+    quiet = TRUE
+  )
+
+  # Invisibly return the file path
+  invisible(output_gpkg)
+}