210127_artikkel/gpkg/01_piirkonnale_ruudustike_lisamine.R

#' Piirkonna ruudustike leidmine
#'
#' Etteantud piirkonna geomeetrilise piirjoone ('piir') järele leitakse selle piirkonna EPSG:3301 projektsioonile vastav piirikast ('boundarybox_3301'). Piirikasti järele leitakse selle kastiga kaetud 5x5 km ('epk10t_grid'), 1x1 km ('epk2t_grid') ja 100x100 meetrit ('epk02t_grid') ruudustikud. Samuti leitakse piirkonna piiri sees olevad ja statistilist ning muud infot sisaldavad 5x5 km ('epk10t') ja 1x1 km ('epk2t') sisaldavad ruudu. Andmed salvestatakse GPKG faili kihtidena. Faili nimeks on objekti nimi.
#'
#' @param obj str Objekti nimi.
#' @param pk Piirkonna geomeetriline joon.
#' @param gpkg_home path Salvestatavate GPKG faili asukoht.
#'
#'
piirkonnale_ruudustike_lisamine <- function(obj = NULL, pk = NULL, gpkg_home = "/tmp") {
  if (is.null(pk) || !sf::st_is_valid(pk)) {
    cat("\nPalun kontrolli geomeetrilise kujundi õigsust.\n")
    return()
  }
  if (is.null(obj)) {
    cat("\nPuudu on objekti nimi.\n")
    return()
  }
  cat(sprintf("\nAlgparameetrid: objekti nimi %s ja GPKG faili kataloog %s\n", obj, gpkg_home))
  dsn <- sprintf("%s/%s.gpkg", gpkg_home, obj)
  input_layer_name <- "piir"
  input_layer <- sprintf("%s|layername=%s", dsn, input_layer_name)
  output_layer_name <- "boundarybox_3301"
  tmp_gpkg_file <- tempfile(fileext = ".gpkg")
  # write to gpkg
  sf::write_sf(pk,
    dsn = dsn,
    layer = "piir", driver = "gpkg", append = FALSE
  )

  ## ------------ 2. piirkonna 3301 projektsiooniga piirikast --------------
  ## 2.1  Esmalt leiame milliset 500x500 meetri ruutudega on objekt kaetud.
  # ruut::qgis_algorithm_search_by_word(str = "grid")
  algorithm <- "native:creategrid"
  # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
  result <- qgisprocess::qgis_run_algorithm(
    algorithm = algorithm,
    CRS = "EPSG:3301",
    # EXTENT = '25.454305528,26.259893095,59.454118579,59.714621582 [EPSG:4326]',
    EXTENT = input_layer,
    HOVERLAY = 0,
    HSPACING = 500,
    TYPE = 2,
    VOVERLAY = 0,
    VSPACING = 500,
    OUTPUT = qgisprocess::qgis_tmp_file(ext = ".gpkg")
    # .quiet = TRUE
  )
  result
  pk_grid <- sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT"))
  # sf::st_geometry(pk_grid) %>% plot()
  ## 2.2 Leitud ruutudele leiame piirikasti.
  algorithm <- "qgis:minimumboundinggeometry"
  # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
  result <- qgisprocess::qgis_run_algorithm(
    algorithm = algorithm,
    FIELD = "",
    INPUT = pk_grid,
    TYPE = 3,
    OUTPUT = tmp_gpkg_file
    # .quiet = TRUE
  )
  result
  pk_grid_bb <- sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT"))
  sf::st_geometry(pk_grid_bb) %>% plot()
  sf::st_geometry(pk) %>% plot(add = T)

  system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s", dsn, tmp_gpkg_file, output_layer_name))

  ## ------------- 3. piirkonna epk10t (5x5 km) ruudud --------------------
  # 3.1 kogu ruutvõrgustik
  ruudud <- c("epk10t_grid", "epk2t_grid", "epk02t_grid")
  j <- 4
  for (j in 1:length(ruudud)) {
    ruut <- ruudud[j]
    output_layer_name <- ruut
    conf <- ruut::get_config()
    pg <- ruut::construct_ogr2ogr_PG_connect_str(conf = conf)
    input <- sprintf(
      'postgres://dbname=\'%s\' host=%s port=%s user=\'%s\' sslmode=%s password=\'%s\' key=\'fid\' srid=3301 type=Polygon checkPrimaryKeyUnicity=\'1\' table=\"%s\".\"%s\" (geometry)',
      conf$dbname, conf$host, conf$port, conf$user, conf$sslmode, conf$password,
      "maaamet", ruut
    )
    # ruut::qgis_algorithm_search_by_word(str = "extract")
    algorithm <- "native:extractbylocation"
    # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
    result <- qgisprocess::qgis_run_algorithm(
      algorithm = algorithm,
      INPUT = input,
      INTERSECT = sprintf("%s|layername=%s", dsn, "boundarybox_3301"),
      OUTPUT = tmp_gpkg_file,
      PREDICATE = c(0, 1)
      # .quiet = TRUE
    )
    result
    assign(ruut, sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT")))
    system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s", dsn, tmp_gpkg_file, output_layer_name))
  }
  sf::st_geometry(epk10t_grid) %>% plot(add = T)

  # 3.2 ainult piirkonnaga seotud ning informatsiooni sisaldav ruutvõrgustik
  ruudud <- c("epk10t", "epk2t")
  j <- 4
  for (j in 1:length(ruudud)) {
    ruut <- ruudud[j]
    output_layer_name <- ruut
    conf <- ruut::get_config()
    pg <- ruut::construct_ogr2ogr_PG_connect_str(conf = conf)
    input <- sprintf(
      'postgres://dbname=\'%s\' host=%s port=%s user=\'%s\' sslmode=%s password=\'%s\' key=\'fid\' srid=3301 type=Polygon checkPrimaryKeyUnicity=\'1\' table=\"%s\".\"%s\" (geometry)',
      conf$dbname, conf$host, conf$port, conf$user, conf$sslmode, conf$password,
      "maaamet", ruut
    )
    # ruut::qgis_algorithm_search_by_word(str = "extract")
    algorithm <- "native:extractbylocation"
    # cat(qgisprocess::qgis_show_help(algorithm = algorithm))
    result <- qgisprocess::qgis_run_algorithm(
      algorithm = algorithm,
      INPUT = input,
      INTERSECT = input_layer,
      OUTPUT = tmp_gpkg_file,
      PREDICATE = c(0, 1)
      # .quiet = TRUE
    )
    result
    assign(ruut, sf::read_sf(qgisprocess::qgis_output(result, "OUTPUT")))
    system(sprintf("ogr2ogr -f GPKG -overwrite  %s %s -nln %s", dsn, tmp_gpkg_file, output_layer_name))
  }
  # sf::st_geometry(epk10t) %>% plot(add = T)
  ## ---------------------- vaata layer'id ----------------------
  # Vaata layer'eid
  sf::st_layers(dsn = dsn)
}
# piirkonnale_ruudustike_lisamine(obj = NULL, pk = NULL, gpkg_home = "/tmp")